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each  region  into  simulation  cells.  Each  cell  is  filled  with  a  number  of 
simulated  molecules  relative  to  the  cell  volume  and  local  number  density. 

The  simulation  is  performed  for  each  region  separately  and  contains: 

*  molecular  motions 

*  generation  of  new  molecules  to  simulate  input  flows 

*  deactivation  of  molecules  to  simulate  output  flows 

*  molecular  collisions 
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used  as  input  data  for  consecutive  runs. 
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ABSTRACT 

Some   procedures    have    been   developed    to   analyze    the   flowfield   of    highly 
underexpanded   axisymmetric   ring  jets    operated   at   high  altitudes.      The  Method 
of    Characteristics    (HOC)    was    used    to   compute    the   Prandtl-Meyer   expansion   fan 
and    the   flow  parameters    in   that   region.      The  MOC  may   also   be   used    to   obtain 
some   indications    about    the   repetitive   expansion   —   compression   behavior   of    the 
jet   as   well   as    the   divergent   shape   of    the   expansion   part   downstream,    when 
the   ambient   pressure   goes    below  certain   limits. 

The   continuum  methods   may   be   used    as    far   as    the   limit  where    trans la tional 
equilibrium  ceases    to   exist.      The    breakdown   of    the   continuum   theory  may   be 
evaluated   using    the   experimental    breakdown   parameter   as    proposed   by 
G.   A.    Bird.      Beyond    this    limit,    the  molecular   Direct-Simulation-Monte-Carlo 
method    (DSMC)    is    applied. 

The   axisymmetric   Monte   Carlo    Simulation    Degins    outside    the   nozzle,    where 
the   breakdown   parameter    has    a  value   of    approximately   0.05.      The   actual   shape 
of    this    breakdown   limit    is    a   closed    lobe   surface   which  starts    at    the   nozzle 
lips,    approximately   follows    a  specific   Mach  wave    in    the    Prandtl-Meyer   fan   and 
closes    back    to    the    axis    far   downstream.      Because    our    interest    is    limited    to 
the   close   vicinity   of    the    corners,    there    the   shape   of    this    surface   may   be 
approximated    by   a  straight    line    (making   a   cone    in   an   axisymmetric    flow). 

For    the   simulation   purposes ,    the   domain   between    the    breakdown   surface   and 
the   wall    is    divided    into   sectors,    each  sector   divided    into    radial    regions    and 
each   region    into   simulation   cells.      Each   cell    is    filled   with   a  number   of 
simulated   molecules    relative    to    the    cell   volume   and    local   number   density. 
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The  simulation  is  performed  for  each  region  separately  and  contains: 

*  molecular  motions 

*  generation  of  new  molecules  to  simulate  input  flows 

*  deactivation  of  molecules  to  simulate  output  flows 

*  molecular  collisions 

Because  there  is  no  apriori  information  about  the  flow  interaction  between 
different  regions,  the  whole  simulation  is  done  on  an  iterative  basis.   A 
first  run  will  supply  the  output  flux  from  each  region.   These  results  are 
used  as  input  data  for  consecutive  runs. 

The  program  runs  as  far  as  the  collision  frequency  is  still  high  and  the 
mean  free  path  is  low  compared  with  the  size  ot  the  cells.   Beyond  this  limit 
the  flow  may  be  regarded  as  collisionless ,  and  the  flux  towards  the  solid  wall 
may  be  computed  directly. 
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I.   INTRODUCTION 

Gas  jets  released  from  spacecrafts  and  external  flows  about  vehicles  at 
high  altitudes  have  a  renewed  interest  in  particular  with  regard  to  two  main 
aspects : 

a.  contamination  of  the  spacecraft  walls 

b.  optical  disturbances  caused  by  the  plume. 

A  spacecraft  gas  dynamics  laser  releasing  a  large  quantity  of  gas  is  highly 
affected  by  these  two  factors  and  the  interest  in  analyzing  them  and  being 
able  to  control  them  have  a  unique  importance  in  further  development. 

A  spacecraft  laser  is  assumed  to  have  a  long  cylindrical  shape  with  the 
optical  power  output  devices  installed  at  one  of  its  bases  as  shown  in  Figure 
1. 

The  output  gas  is  released  through  a  ring  nozzle,  undergoes  a  fast 
three  dimensional  axisymmetric  expansion,  and  forms  a  plume  covering  the  whole 
meridian  plane  of  the  vehicle.   It  widens  to  large  angles  of  expansion  so  that 
it  may  intersect  the  laser  beam.   Back  scattered  molecules  may  return  to  the 
wall  of  the  spacecraft  causing  contamination  and  degradation  of  surfaces  and 
vehicle  parts. 
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Figure    1.      The    Spacecraft    Laser. 

M     -  Mach   number    at    the   exit    surface      ~  4     the    jet    gas    is    composed   of    two 

species 

heavy  molecules      Mqj   =    19 

light    molecules      I%2  =      ^ 
Altitudes    between   200  to    1000   km. 


Continuum  flow  theory  may  be  used  to  solve  the  flowfield  and  flow 
parameters  as  far  as  there  is  translational  equilibrium,  it  means  that 
intermolecular  interactions  are  fast  enough  to  maintain  expansion  rates. 
Wherever  these  interactions  are  too  slow,  the  continuum  flow  becomes  invalid 
and  the  molecular  flow  theory  should  be  employed. 

The  solution  for  the  continuum  regime  is  computed  here  by  means  of  the 
Method  of  Characteristics  (HOC)  [1,2,3].   The  limit  where  continuum  breakdown 
occurs  was  estimated  by  the  experimental  breakdown  parameter  as  proposed  by  G. 
A.  Bird  [4].   Beyond  this  limit,  it  is  proposed  to  compute  the  molecular  flow 
by  means  of  the  Direct  Simulation  Monte  Carlo  technique  as  described  in  detail 
by  Bird  [4] . 

For  moderate  and  low  pressure  ratios  (static  pressure  at  the  nozzle  exit 
to  the  ambient  pressure),  an  underexpanded  jet  exhibits  a  repetitive  expansion 
-  compression  behavior  with  a  geometry  depending  on  the  initial  Mach  angle,  on 
the  Prandtl-Meyer  fan  angle,  and  on  the  gas  specific  heats  ratio.   For  lower 
amDient  pressure  which  occurs  at  higher  altitudes,  the  first  compression 
region  is  pushed  out  to  the  envelope  of  the  jet  forming  the  barrel  shock.   If 
the  ambient  pressure  is  low  enough,  this  compression  region  may  disappear  due 
to  the  molecular  behavior  of  the  flow  [6,11]. 

The  breakdown  of  the  continuum  theory  occurs  in  a  region  where  the  gas 
density  and  pressure  are  high  compared  with  the  density  and  pressure  in  the 
ambient  gas.   At  high  altitudes  the  ratios  between  these  parameters  may  reach 
10°  or  more.   Considering  this  range  of  variations,  computational 
validity  dictates  the  use  of  the  Direct  Simulation  Monte  Carlo  method.   In  the 
higher  density  range  the  jet  will  be  considered  as  consisting  of  two  species 
of  gas,  their  molecular  model  will  be  "the  hard  sphere  molecule"  model  and 


ambient  gas  is  not  allowed  to  protrude.   In  the  lower  density  region  the  flow 
will  be  regarded  as  collisionless . 

In  the  following  chapters  we  bring  the  detailed  description  of  the 

computer  programs  which  solve  the  different  parts  of  the  flowfield. 


II.   THE  CONTINUUM  REGIME 

A.   THE  TWO  DIMENSIONAL  ISENTkOPIC  UNDEREXPANDED  JET 

The  results  brought  here  are  based  on  the  supersonic  steady  isentropic 
flow  theory  as  described  in  literature  (see  for  example,  Shapiro  [1],  Liepman 
and  Roshko  [2],  and  Owczarek  [3]). 

The  ranges  of  paremeters  of  a  jet  emerging  from  a  gas  dynamics 
spacecraft  laser  are: 

a.  The  Mach  number  at  the  exit  surface  Mq  «  4.   The  static  pressure 
at  the  exit  plane  P0  »  136pa.   Ambient  pressure  (Pamb)>  temperature 
(Tamb)  and  other  thermodynamic  properties  of  amoient  gas  depend  on 
the  altitude  as  shown  in  Table  1.   The  jet  gas  may  consist  of  DF ,  HF , 
Helium  and  other  species.   In  the  programs  we  limit  the  compositon  to 
two  species:   Air  and  He.   (The  program  allows  changes  in  the 
composition  and  types  of  gas.) 

P 

b.  The  pressure  ratio  — —  for  the  minimum  required  altitude  (200 

amb 

km)  assures  that  the  jet  is  highly  underexpanded  (we  show  later  the 

influence  of  this  ratio  on  the  shape  of  the  jet). 

The  following  thermodynamic  relations  are  valid  as  long  as  the 

compessible  flow  is  isentropic 

TT  =  T(l  +^~  M2)  (1) 

JL 


PT  =  P(l  +  JCl  m2)^1  U) 
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,PT  -.p(l  -  JCi  M2)^1  (3) 
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where   T-p,  Pt»  an^  PT  are  tne  total  temperature,  pressure  and  density 
(constant  for  isentropic  field).   T,  P,  and  p   are  local  temperature,  pressure 
and  density  y     is  the  specific  heat  ratio  of  the  gas  (considered  here  as 
constant),  11  is  the  local  Hach  number. 

The  partial  differential  equation  of  motion  for  supersonic  2-D 
irrotational  and  isentropic  flow  is  a  hyperbolic  equation  having  solutions 
obtained  from  invariants  along  characteristic  lines.   Physical  interpretation 
of  these  lines  are  the  compression  or  expansion  waves  which  are  oriented  at 
Mach  angles  relative  to  the  streamlines. 

Once  the  directions  of  the  characteristics  (, waves)  are  determined 
everywhere  in  the  field,  all  other  parameters  may  be  calculated. 

B.   THE  TWO  DIMENSIONAL  PLANAR  JET 

Tne  compressible  supersonic  jet  flow  is  characterized  by  two  families 
of  characteristics  (pressure  waves)  starting  at  each  corner  of  the  nozzle 
lips.   Each  of  these  families  of  waves  forms  a  Prandtl-lieyer  fan.   The 
streamlines  crossing  the  characteristic  waves  bend  outwards  resulting  in  an 
increase  in  the  flow  area.   The  angle   u   between  the  streamline  and  the 
pressure  wave  is  a  function  of  the  local  ilach  number  as 

u  =  arcsin  (1/M)  (4) 

The  symbol   u   is  the  Ilach  angle. 


Using    the    isentropic   relations,    we   can  find    the   relation   between    the 

turning   angle    ( 9)    and    the   local   Mach  number    (M)    as 

2 
(M     -    1) 

d9  =  -     *= r?~  dM  (5) 

M(  1   +  ^j-  IC ) 

Integration   between    the   conditions      M=l      and   M     gives    the    total    turning 
angle   starting   at    the    throat    (where     M=l)    up    to   a  point  with  given     M.      The 
result  gives    the   Prandtl-Meyer   function    (angle)    as 

v(M)    =     ^-       arctg     ^-  (M2-l)      -     arctg    (M2-l)  (6) 

In  the  close  vicinity  of  the  nozzle  lips  where  the  two  families  of 
characteristics  do  not  intersect  with  each  other  there  is  a  "simple  region"  of 
expansion.   There  the  flow  parameters  are  defined  by  v  and  9  of  each 
characteristic  line.   Further  downstream  the  waves  intersect  each  other.   In 
this  part  of  the  flow  parameters  are  defined  by  the  two  intersecting 
characteristics.   A  singularity  occurs  when  the  initial  Mach  number  is  unity 
and  a  special  treatment  is  required  to  start  the  calculations  at  that  point 
(this  special  treatment  has  not  been  brought  here).   A  particular  importance 
of  the  Prandtl-Meyer  function  is  when  analyzing  the  two  dimensional  flow  using 
the  hodograph  plane. 

C.   THE  HODOGRAPH  PLANE  FOR  A  TWO  DIMENSIONAL  JET 

The  hodograph  plane  is  a  representation  of  the  flow  parameters  in  the 
velocity  plane.   Figure  2  shows  a  hodograph  plane  calculated  for  a  simple 
gas  (y=1.4).   The  circles  represent  constant  Mach  numbers  and  constant 


a*      v  ■- — 


<r-/ 


9]^",    9i+  -  the    initial   turning 


angles  for  Mn>l 
0^  ,  62  ~  maximum  turning 
angles  for  a  highly  under- 
expanded   jet 

8_,v  -  maximum  turning   angle 
for    the   specific    gas. 


Figure    2.      The   Hodograph   Plane. 


1-2-3-4-1    defines    a  jet   with     MQ>1   expanding   into      PamDl   <   p0» 
resulting  simple   repetitive   expansion-compression. 
l*-2*-3*-3**-4*-l*  defines    a  jet    with  Mq-1      expanding   into 
Pamb2   ^   Po      *    resulting  in   a  highly   underexpanded   jet. 


p 
pressure  ratios   (-5—)   •   The  epicycloids  are  the  two  families  of 

T 
characteristics.   The  angles   9  are  the  turning  angles  due  to  expansion  or 
compression  (which  are  both  present  in  supersonic  jets). 

To  define  an  isentropic  supersonic  jet  on  the  hodograph  plane  it  is 

necessary  to  define  the  Mach  number  at  the  exit  plane,  the  pressure  ratios 

p 
p            cirnb 
(— — )   and   (— )  ,  and  the  specific  heat  ratio  ( y)  of  the  gas. 

T  T 

Investigating  the  shapes  of  the  jet  as  a  function  of  ranges  of  parameters 
we  may  get: 

Pamb 

a.  simple  underexpanded  jets  for  which  — is  high  enough  so  that 

T 
the  two  families  of  characteristics  intersect  each  other. 

b.  a  critical  underexpanded  jet  for  which  — is  low  enough  so  that 

T 

the  intersection  between  the  outer  characteristic  lines  of  the  two 
familes  occur  on  the  outer  hodograph  circle. 

c.  highly  underexpanded  jets  for  which  a  part  of  characteristics  do  not 
intersect  at  all. 

d.  expansion  into  complete  vacuum  so  that  there  are  no  reflections  from 
the  jet  boundaries  and  therefore  the  compression  region  disappears. 
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D.   THE  SHAPES  OF  TWO  DIMENSIONAL  JETS 

The  following  paragraphs  further  detail  the  different  shapes. 
1 .   Simple  Underexpanded  Jets 

Figure  3  shows  the  physical  plane  and  the  hodograph  plane  of  a  simple 
underexpanded  jet.   If  IIqM  ,  an  initial  turn  of  the  flow  90  is  made  within 
the  nozzle.   An  aditional  turn  of  8  is  due  to  the  underexpansion.   9  is  found 
by  the  intersection  of  characteristics  (1-2)  with  the  circle  defined  by 

Pamb/PT  • 

When  M0>1  ,  the  characteristic  line  1-2  is  described  in  the  physical 
plane  by  a  region  in  which  only  one  family  of  expansion  waves  are  present 
(simple  region,  see  transverse  line  between  1  to  2  in  the  physical  plane).   A 
different  family  of  expansion  waves  forms  a  second  simple  region  when  moving 
between  2  to  3  .   At  a  larger  distance  from  the  exit  plane,  reflected  waves 
from  the  free  streamline  (jet  boundary)  cause  compression.   An  ideal 
representation  of  such  a  jet  is  a  repetitive  pattern  of  expansion  and 
compression. 

For   M0=l  ,  the  tail  waves  of  both  families  are  perpendicular  to 
the  flow,  thus,  both  lie  on  line  AB  .   That  means  that  if   M0=l   there  is  no 
simple  region  near  the  exit  plane  of  the  jet.   The  characteristic  line  1-2  on 
the  hodograph  plane  becomes  a  single  point   A  (or  B)  when  located  in  the 
physical  plane. 

2  .   Critical  Underexpanded  Jets 

We  define  a  critical  underexpanded  jet  when  point  (3)  (see  Figure  4) 
lies  on  the  limiting  circle  M^,  or  P/P^O.   This  means  that  there  is  a  core 
within  the  jet  where  the  pressure  approaches  zero  and  Uach  approaches 
infinity.   This  core  is  theoretically  bounded  at  its  upstream  side  by 
expansion  waves  and  downstream  by  compression  waves. 
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Jet  boundary  (freestream) ,  P  =  Pamh 
L, 


compression  waves 


. ,  expansion  waves 


streaml ine 


©  uniform  axial    flow   (P=PQ) 

(D  uniform  flow  parallel    to   Lo(p=pamb^ 

<D  uniform  axial    flow   (lowest   pressure-Pmin 

@  uniform   flow  parallel    to   l-4(p=pamD) 

(D  uniform  axial    flow   (P=PQ) 


Figure    3.      Flow  at    exit    of    a  simple    under  expanded   jet 

(a)  physical  Plane 

(b)  Hodograph  Plane 
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Using  the  theoretical  expressions  for  isentropic  ideal  flow  one 

P       P 

may  derive  the  values  of  — or  — as  function  of  M   which  causes  a  jet 

Pt      Pq 

to  be  critical  underexpanded. 

P  amb 

Figure  5  shows  results  of   as  functions  of  M0  for  gases   with 

o 

different  specific  heat  ratio. 

3.   Highly  Underexpanded  Jets 

When  Pamb   ^s  lower  than  the  critical  values  as  shown  in  figure  4, 
point  3  does  not  exist  (there  is  no  intersection  between  lines  2-3  2'-3'). 
This  means  that  the  repetitive  reversible  expansion/compression  shapes  ceases 
to  exist.   The  envelope  of  the  jet  starts  at  an  angle  defined  by   9  at  the 
nozzle  exit  plane  and  approaches  an  asymptotic  angle  defined  by    Q±im      (see 
Figure  6). 

In  this  case  we  get  two  (symmetric)  groups  of  characteristics  C^-C2  and 
C'^-C'2  (Figure  6)  with  no  intersection  between  them.  C±    and  Cj '  define  the 
inner  limit  for  reflected  characteristics,  C2  and  C2 '  define  the  outer  limit. 
As  the  rest  of  reflected  waves  lie  between  C2  or  C'2  and  the  jet  boundary,  we 
may  conclude  that  a  compression  region  may  exist  only  in  a  layer  along  the  jet 
boundary. 

because  of  irreversible  effects  such  as  shear  stresses,  heat  transfer  due 
to  high  temperature  gradients,  or  condensation  effects  in  real  gases,  the 
compression  layer  may  be  interpreted  as  the  "barrel  shock". 
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Figure  4.   The  Hodograph  Plane  for  a  Critical  Underexpanded  Jet. 
(1-2-3-4-1),  (t-1.4) 
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Figure    5.      Dependence    of   critical   values    of      pamb/po 

on  the   Mach   number    at    the   exit    plane      (MQ)      for    different    values 
of    specific   heat    ratio    (y). 
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Figure  6.   The  Hodograph  Plane  for  a  Highly  Underexpanded  Jet 


(Y  -    1.4) 

9max  -  maximum  turning  angle 

6         -  total   turning   angle    in  the    specific   jet      9  > 

Q±±m  -  limiting   angle    of   compression  region. 


3        /o 
'max'  2 
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Figure  (7)  shows  the  hodograph  plane  for  an  air  jet  (^=1.4)  expanding 
into  an  ambient  pressure  105  times  lower  than  the  static  pressure  at  the  exit 
plane.   Figure  (8)  is  a  schematic  description  of  the  jet  (the  data  for 
this  jet  is  given  in  Table  2).   The  shape  shown  in  Figure  (8)  may  be  compared 
with  the  barrel  shock  photograph  in  page  208  Reference  [3]. 

4.   EXPANSION  INTO  A  COMPLETE  VACUUM  -  Pamb=Q 

This  is  an  extreme  situation  in  which  the  maximum  turn  angle  of 
streamlines  occur  near  the  nozzle  exit  plane.   The  theoretical  free  stream  is 
defined  only  by  the  theoretical  turning  angles.   All  streamlines  in  the  flow 
expanded  monotonically  towards   P«0  without  being  reflected  by  the  jet 
boundaries . 
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TABLE  2 


CALCULATION  OF  CHARACTERISTICS  ON  PHYSICAL  PLANE 


State 

C 

9* 

8 

v(M) 

M 

U 

9+u 

9~U 

P/Pt 

1 

0 

0 

0 

0 

1 

90° 

90° 

-90° 

0.5283 

2 

0 

10 

10 

10 

1.435 

44.15 

54.15 

-34.15 

3 

0 

20 

20 

20 

1.775 

34.3 

54.3 

-14.3 

.2990 

4 

0 

40 

30 

30 

2.135 

27.93 

57.93 

2.07 

.1810 

5 

0 

60 

40 

40 

2.54 

23.18 

63.18 

16.82 

.1035 

6 

0 

80 

50 

50 

3.013 

19.38 

69.38 

30.62 

.55   -1 

7 

0 

100 

60 

60 

3.595 

16.15 

76.15 

43.85 

.2672  -1 

8 

0 

120 

67.7 

67.7 

4.145 

13.96 

88.7 

53.74 

.1146  -1 

9 

20 

135.4 

0 

20 

1.775 

34.3 

34.3 

-34.3 

.5437  -2 

10 

20 

20 

10 

30 

2.135 

27.93 

37.93 

-17.93 

.1810 

11 

20 

40 

20 

40 

2.54 

23.18 

43.18 

-  3.18 

.135 

12 

20 

60 

30 

50 

3.013 

19.38 

49.38 

10.62 

.55   -1 

13 

20 

80 

40 

60 

3.595 

16.15 

56.15 

23.85 

,2672  -1 

14 

20 

100 

50 

70 

4.339 

13.315 

63.32 

36.69 

.1146  -1 

15 

20 

120 

57.7 

77.7 

5.085 

11.34 

69.04 

46.4 

.425   -2 

16 

40 

135.4 

0 

40 

2.54 

23.18 

23.18 

-23.18 

17 

40 

40 

10 

50 

3.013 

19.38 

29.38 

-  9.38 

.55   -1 

18 

40 

60 

20 

60 

3.595 

16.15 

36.15 

3.85 

.2672  -1 

19 

40 

80 

30 

70 

4.339 

13.315 

43.32 

16.69 

.1146  -1 

20 

40 

120 

40 

80 

5.3479 

10.777 

50.78 

29.23 

.4215  -2 

21 

40 

135.4° 

47.7 

87.7 

6.433 

8.943 

56.64 

38.76 

22 

60 

60 

0 

60 

3.595 

16.15 

16.15 

-16.15 

.1146  -1 

23 

60 

80 

10 

70 

4.339 

13.315 

23.32 

-  3.32 

.4215  -2 

24 

60 

100 

20 

80 

5.3479 

10.777 

30.78 

9.22 

25 

60 

120 

30 

90 

6.8190 

8.433 

38.43 

21.57 

26 

60 

135.4° 

37.7 

97.7 

9.2105 

6.2330 

43.93 

31.5 

27 

80 

80 

0 

80 

5.3479 

10.777 

10.777 

-10.77 

28 

80 

100 

10 

90 

6.8190 

8.433 

18.433 

1.57 

29 

80 

120 

20 

100 

9.2105 

6.233 

26.23 

13.8 

30 

80 

135.4 

27.7 

107.7 

12.45 

4.61 

32.3 

23.1 

31 

100 

100 

0 

100 

9.2105 

6.233 

+  6.233 

-  6.233 

32 

100 

120 

10 

110 

13.874 

4.1331 

14.13 

5.9 

33 

100 

135.4 

17.7 

117.7 

21.4 

2.678 

20.37 

15.0 

34 

120 

120 

0 

120 

27.335 

2  097 

2.097 

-  2.097 

35 

120 

135.4 

7.7 

127.7 

104 

0.546 

8.2 

7.2 

36 

- 

- 

11° 

~ 

00 

0 

- 

- 

0 

51 

20 

115.4 

47.7 

67.7° 

4.145 

13.96 

33.7 

.5437  -2 

52 

40 

115.4 

37.7 

77.7 

5.085 

11.34 

26.36 

53 

40 

95.4 

+27.7 

67.7° 

4.145 

13.96 

13.7 

.5437  -2 

54 

60 

75.4 

+  7.7 

67.7° 

4.145 

13.96 

-  6.26 

.5437  -2 

55 

60 

95.4 

17.7 

77.7 

5.085 

11.34 

6.4 

56 

60 

115.4 

27.7 

87.7 

6.433 

8.943 

18.8 
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Figure  7.   Hodograph  plane  showing  Che  points  on  Che  mesh  of  characteriscics 
(RelaCed  Co  physical  plan  Figure  8). 
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Figure  8.   Schematic  shape  of  a  highly  underexpanded  jet 
(See  Figure  7  -  The  Hodograph  Plane) 
The  arrows  indicate  flow  direction. 
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E.      THE   METHOD   OF   CHARACTERISTICS;    COMPUTATION   OF   PLANAR  AND   AXISYMMETRIC 
TWO-DIMENSIONAL   FLOWS 

For   a  planar   two  dimensional  flow,    the   Prandtl -Meyer   function   v  and 
flow  direction    9  at    any  point    (3)   in  the   field  may   be   calculated  using  data  of 
two  other   points    (1)   and    (2)   located  on  characteristic   lines   that    intersect    at 
(3)    (see   Fig.    (9)) 


1 


1 


v3  ■  "J  (v1  +   v2)   +  ~2   (81   -    92) 


e3  --j  (vx  -  v2)  +  ±  (9l  +  e2) 


(7) 

(8) 


e-  y2 


Figure    9.      The   calculation   of    v  and    9  for   point    (3)    is    based   on   data   for 
points    1   and   2. 
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For   axjsymmetrlc ,   two-dimensional  flow,    Liepmann  and  Roshko    [2]    developed 
expressions   for   finding  the   Prandtl-Meyer   function   (v)   and  the   flow  direction 
(8)  which  are   given  by: 


sine  sine, 

v3  -1  <VL  +   v2)   +  \  (9l   -   e2)   +I  [sin   ^  _  ACn  +  sin   ^ 


F7~  An23] 


(9) 


1                             1                             1                    Sin9l 
93  =  I  (V1   '   V  +1  (91  +  V   +I  [sin   pl  "7 AC13  "  Sin   y 


sin62 
2"rT"  An23] 


The   angles   and  subscripts    are   shown  in  Figure    (10). 


(10) 


fln/s   of  2>yMna.TRY 


Figure    10.      Calculation   of    6  and    v  for    axisymmetric   flow, 
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It  is  obvious  that  for  the  axially  symmetric  flows  the  increase  in  the 
radius  causes  the  increase  in  the  flow  cross  section  and  influences  the  flow 
direction  and  the  Prandtl-Meyer  function.   These  facts  have  been  taken  into 
account  when  developing  the  "compatibility  equations"  (9,10). 

Using  equations  (9,10)  we  have  developed  a  computer  program  which 
enables  the  calculation  of  the  jet  flow  for  a  two  dimensional  and  for  an 
axially  symmetric  geometry  (ring  jet). 

The  listing  of  the  program,  the  program  description  and  some  results  are 
given  in  Appendix  A. 
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III.   THE  BREAKDOWN  OF  THE  CONTINUUM  THEORY 

A.   GENERAL  CRITERIA 

As  described  in  detail  by  Bird  (Chapter  1  Reference  [4]),  the  validity 
of  the  continuum  approach  has  been  identified  with  the  validity  of  the  Navier 
Stokes  equations.   This  requires  that  the  Knudsen  number  Kn  =  A/L  should  be 
small  compared  with  unity  (X  is  the  mean  free  path  and  L  is  a  scale  length  for 
the  specific  flow  field).   For  Kn  larger  than  a  certain  limit  (between  0.01 
to  0.1  depending  on  the  required  accuracy)  a  microscopic  approach  is 
necessary. 

For  small  values  of   L  ,  the  microscopic  approach  may  lead  to  statistical 
fluctuations  of  the  results  due  to  the  small  number  of  molecules  participating 
in  the  flow  processes.   In  figure  (11)  which  was  reproduced  from  Bird's  book 
[4,  figure  1.6],  the  regimes  of  rarefied  flow  and  high  fluctuations  are 
depicted.   The  flow  around  the  jet-air  boundaries  near  a  spacecraft  is 
generally  rarefied  (high  Knudsen  number),  but  has  insignificant  fluctuations. 
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Figure    11.      Limits    for    continuum  approach   and   microscopic    approach 
(d-3.7   x   10-10  m). 
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B.   THE  EMPIRICAL  CRITERION 

The  Method  of  Characteristics  (MOC)  was  used  to  compute  the  jet  flow 
and  the  results  obtained  from  the  computer  program  "AXSYM"  are  valid  as  long 
as  the  continuum  flow  theory  is  valid. 

The  continuum  flow  requires  that  the  mean  free  path  should  be  negligibly 
small  in  comparison  with  the  scale  length  of  the  macroscopic  flow  variations. 
The  classical  theory  for  Prandtl-Meyer  expansion  may  therefore  be  expected  to 
fail  at  progressively  larger  distances  from  the  nozzle  lips  as  the  gas  density 
decreases  with  the  increasing  flow  angle  and  Mach  number.  The  empirical 
criterion  for  the  breakdown  of  continuum  flow  in  steady  expansion  flow  [4 J  is 
that 


P     S-* 


pv 


dp_ 

dS 


=  0.05 


(11) 


where 


q   =  stream  velocity 

p    =  density 

v  =  molecular  collision  frequency 


dfi.1  _ 
dS  I 


absolute  change  in  density  while  moving  a  distance 


dS   along  a  streamline 
Introducing  the  breakdown  parameter  P  into  Che  program,  gives  the 
definition  of  the  boundary  where  the  flow  should  be  calculated  Dy  means  of  the 
molecular  flow  theory,  i.e.,  by  solving  the  Boltzmann  equation. 
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For  an  underexpanded  jet  with  a  high  initial  Mach  number, 
the  breakdown  surface  is  nearly  a  streamline.   Furthermore,  the  range  of  flow 
parameters  for  the  present  problem  are  such  that  the  simple  region  extends  to 
very  large  distances  and  near  the  nozzle  lip  the  breakdown  limit  may  be 
approximated  by  a  straight  line. 

For  the  axisymmetric  jet  there  is  no  simple  region,  however,  for  the 
region  of  interest  it  may  be  regarded  as  linear. 

The  method  proposed  in  the  present  work  for  solving  the  flow  behind  the 
breakdown  boundary  is  the  Direct  Simulation  Monte  Carlo  (DSMC).   For  tnis 
purpose,  a  computer  program  "SIMUL"  was  developed.   In  the  following  chapter 
we  describe  the  algorithms  required  for  the  specific  problem,  the  geometry  and 
the  data  organization.   Detailed  program  description  is  given  in  Appendix  (B) . 
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IV.   THE  MOLECULAR  FLOW  IN  AN  AXISYMMETRIC  RING  JET 

A.   GENERAL  CONSIDERATIONS 

The  part  of  the  field  in  which  the  jet  may  be  calculated  by  means  of  the 
continuum  theory  was  described  in  Chapter  II.   There  we  calculated  also  the 
boundaries  where  continuum  theory  becomes  invalid  and  molecular  calculation 
should  be  employed.   In  fact,  the  molecular  theory  and  the  molecular  Boltzmann 
equations  are  universal  and  hold  for  the  entire  flowfield.   However, 
computational  requirements  make  the  Boltzmann  equation  impractical  for  the 
upstream  flows.   Therefore  we  limit  our  solution  only  to  the  part  of  the  flow 
beyond  the  region  where  continuum  breakdown  occurs. 

As  a  result  obtained  from  MOC  solution  the  "breakdown",  i.e.,  the  locus 
where  the  breakdown  parameter  p   has  values  between  0.03  to  0.06,  for  the 
region  close  to  the  nozzle  lips  this  boundary  may  be  approximated  by  a 
straight  line  (for  axisymmetric  flow  this  line  is  the  envelope  of  a  cone,  see 
Figure  (12)). 

For  the  specific  jet  and  gas,  the  breakdown  occurs  in  a  region  where  the 
number  density  is  in  a  range  of  102*  molecules /m3 .   For  ambient  gas  at  an 
altitude  of  200  km  the  number  density  is  10  ^  and  decreases  to  a  range  of  10  *  1 
at  1000  km.   In  order  to  be  able  to  express  this  vast  change  in  a  simulation, 
we  would  need  to  have  an  extremely  large  number  of  cells  which  would  be 
impossible  to  store  in  a  computer.   To  overcome  this  problem  we  are  required 
to  make  a  less  exact  formulation  which  enables  the  production  of  results, 
having  to  pay  the  penalty  of  "smearing"  the  steep  gradients  and  obtaining 
averages  within  layers  of  simulated  cells.   Unfortunately  it  is  impossible  Co 
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Figure  12.   Regions  in  a  ring  jet. 
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predict  how  far  the  simulated  results  will  be  from  the  exact  solutions.   These 
comparisons  have  to  be  made  after  getting  final  results  of  this  simulation. 
1 .   The  Direct  Simulation  Monte  Carlo  Method 

The  direct  simulation  Monte  Carlo  Method  is  a  technique  for  a  computer 
modeling  of  a  real  gas  by  some  thousands  of  simulated  molecules.   The  velocity 
components  and  the  coordinates  of  the  simulated  molecules  are  stored  in  the 
computer  and  are  modified  with  time  as  a  result  of  collisions  and  boundary 
interactions.   A  detailed  description  of  some  problems  and  their  solutions 
by  means  of  direct  simulation  is  given  in  [4], 

To  follow  the  molecular  motion  it  is  necessary  to  divide  the  simulated 
domain  into  a  network,  of  cells.   The  size  of  a  cell  mus  t  be  such  that  the 
cnange  in  flow  properties  across  each  cell  is  small.   The  time  is  advanced  in 
discrete  steps  DTM,  such  that  DTM  is  small  compared  with  the  mean  collision 
time  per  molecule.   If  there  is  a  flow  going  through  the  domain,  DTM  should  De 
small  compared  with  the  mean  time  required  for  the  mean  flow  to  cross  trie 
cells.*   Both  cell  size  ( ( DR) , ( DDALFA)  -  radial  size  and  angular  size  as  they 
appear  in  the  program)  and  DTM  may  vary  in  the  simulation  with  position  and 
time. 

Applications  such  as  free  jet  expansion  in  which  large  gradients  of  flow 
properties  are  expected,  may  require  a  very  large  number  of  cells  for  the 
simulation.   In  these  cases  the  computer  memory  requirements  to  store  cells' 
data  and  molecules'  data  may  exceed  the  available  computer  storage.   A 


*If  DTM  is  chosen  to  be  very  small  compared  with  the  mean  time  between 
collisions  then  the  simulation  will  require  a  very  large  number  of  runs  such 
that  the  number  of  collisions  will  be  sufficient.   If  DTM  is  large  the 
molecules  are  washed  out  by  ttie  mean  flux  and  there  is  no  time  for  the 
collisions  to  influence  tne  flow. 
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solution  for  this  problem  is  to  divide  the  simulation  space  into  smaller 
regions  and  to  run  the  simulation  for  each  region  separately.   If  there  is  an 
interaction  between  different  regions,  which  apriori  is  undefined,  the 
solution  should  be  found  iteratively.   (That  means  that  each  run  will  provide 
data  for  consecutive  runs  and  the  procedure  should  be  repeated  until  the 
results  converge  to  a  steady  solution. 

The  computation  of  a  representative  set  of  collisions  based  on  mean 
collision  time  per  molecule  is  invalid  for  a  computerized  simulation  because 
of  the  large  computer  time  and  computer  memory  requirements.   Instead,  the 
method  proposed  by  Derzko  which  is  described  in  details  by  Bird  [4]  may  be 
employed.   Following  this  method,  an  averaged  mean  time  between  collisions  of 
species   L   with  species   M  for  a  cell  is  calculated.   The  number  of 
collisions  of  each  type  (L.M  species)  is  such  that  the  collision  time 
counters  are  kept  concurrent  with  the  overall  time  parameter.   The  L.M 
collision  time  for  a  cell  containing   N^  and  Njj  molecules  with  collision 
cross  section   a^i  ,  number  densities  n^   and  nyL     and  relative  velocity   Cr, 
is  given  by 

.,        .LP  1  .MP        _1 ,19x 

At    =  —    —    +  —   ~  ci^; 

c        NT       oT,,n7IC  N..      aTtrnTC 

L        LM  Mr  M        LM  L   r 

where      LP      and   MP      are    the    probabilities    that    the    collision   will    be   effective 

for    the      L      and      M     molecules    respectively. 
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B.   THE  GEOMETRY  OF  THE  SIMULATED  DOMAIN,  SECTORS,  REGIONS  AND  CELLS 

Figure  (13)  shows  a  cross  section  of  the  simulated  domain  for  the 
axisymmetric  (ring)  flow.   Points  A  and  A'  are  the  nozzle  lips.   Starting  at 
"A"  and  assuming  the  "breakdown"  boundary  to  be  a  straight  line,  we  obtain  the 
cross  section  of  the  molecular  domain  as  a  sector  defined  by  LAM.   The  solid 
wall  is  defined  by  AL.   The  arc  LM  may  be  assumed  to  be  far  enough  so  that  the 
pressure  along  it  may  be  assumed  to  equal  the  ambient  pressure.   Molecules 
originated  in  the  jet  cross  the  breakdown  boundary  with  a  velocity,  direction, 
temperature  (and  other  thermodynamic  properties)  as  found  from  the  continuum 
solution. 

The  molecular  domain  LAM   is  divided  into  secondary  sectors,  and  each  of 
these  are  divided  into  several  radial  regions  making  the  "simulation 
regions" . 

Because  we  have  no  apriori  information  on  how  the  expansion  occurs,  the 
angle  of  each  sector  (v/hich  mainly  is  in  the  direction  of  the  expansion 
gradients)  is  left  to  be  a  result  of  the  internal  calculation. 

Each  region  is  divided  into  NRD  radial  divisions  and  NAD  angular 
divisions  making  a  network  of  NAD*NRD  simulation  cells.   The  angle  DALFA  of 
all  regions  in  a  sector  is  constant.   Taking  NAD  constant  for  all  regions  in 
a  sector,  we  get  the  angle  of  a  cell  DDALFA  constant.   Defining  the  raaial 
size  of  a  cell  DR  as  constant  we  get  a  cell  cross  section  area  proportional 
to  the  radius  R  measured  from  the  nozzle  lip  (point  A). 
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Figure    13.     Cross   section  of   the   simulation  domain,    definition  of   sectors, 
regions,    cells    and   coordinates. 
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The  size  of  a  cell:   In  order  to  get  accurate  simulated  results  it  is 
recommended  to  define  the  size  of  a  cell  (DR  and  R*DDALFA)  small  compared  with 
the  mean  free  path  of  the  molecules   A  (typical   DR=A/3).   However,  as  we  do 
not  expect  to  get  large  changes  in  flow  parameters  along  the  radius  we  may 
allow  DR  be  much  larger  than   A/ 3.   The  angular  size  of  the  largest  cell  in  a 
sector  should  comply  with  this  requirement,  but  because  of  the  computer 
limitation  it  is  set  to  be  equal  to  5*  A.   This  will  be  the  basis  for  defining 
DALFA  for  each  sector. 
C.   INITIAL  NUMBER  OF  MOLECULES  IN  CELLS 

A  "reasonable"  number  of  molecules  in  a  simulation  is  several  thousands 
(a  larger  number,  which  is  better,  may  be  used  for  simple  problems  or  when 
using  a  single  user  computer  with  large  user  memory  space).   The  initial 
setting  of  molecules  in  cells  is  usually  based  on  a  guess  of  the  number 
density  in  the  specific  cell.   (The  number  of  molecules  in  cells  will  change 
during  the  simulation  according  to  the  input/output  calculated  fluxes  to  the 
specific  region). 

The  number  density  and  the  size  of  a  cell  are  specified  only  in  three 
dimensional  flows.   When  applied  to  a  two  dimensional  flow  the  simulation  may 
be  regarded  as  applying  to  an  arbitrary  thin  slice  of  the  real  flow.   In  the 
axisymmetric  flow  we  define  the  width  of  a  cell  by  the  angle  DFI  as  shown  in 
Figure  (14),  constant  within  a  region. 

The  initial  number  density  in  cells  of  a  given  region  is  set  constant. 
Defining  the  total  number  of  simulated  molecules  in  the  region  the  number  of 
molecules  in  each  specific  cell  becomes  a  function  of  DFI.   For  example, 
assume  that  we  limit  the  number  of  molecules  in  the  smallest  cell  in  the 
region  to  15  then: 
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Figure    14.     Variation  of   the   angle   DFI. 

To  maintain  the    number    of    simulated  molecules   within  computational 
limits,    the    'width'    of  each  region  defined   by  DFI  is   such  that: 
MIN  ■   number    of   molecules    in  smallest    cell   in  a  region 
MIN  =   VOLUME    (smallest    cell)    *  number    density 

-   f(R,    ALFA,    DALFA)    *  DFI   *  number    density   for    flux 
calculations 
DFI   is    a  weighting   factor. 
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DFI*(contant)  *  (number  density)  =  15 
Other  cells  contain  the  initial  number  of  molecules  proportional  to  their 
volumes . 
D.   DEFINITION  OF  INPUT  AND  OUTPUT  FLOWS  FOR  A  REGION 

The  cross  sections  of  all  regions  (except  those  regions  near  the  nozzle 
lips)  are  quadrilateral.   Through  the  sides  of  the  region  molecules  are 
allowed  to  enter  or  to  leave  according  with  the  boundary  conditions  or  as  a 
result  of  molecular  velocity.   For  the  first  sector,  near  the  breakdown 
boundary  the  input  flow  (FWP1  and  FWP2  see  Figure  (15))  is  defined  by  the 
results  from  the  continuum  flow.   FEN1 ,  FNN1 ,  FSN1 ,  etc.  are  results  of 
counting  and  averaging  the  outgoing  molecules  (the  different  vector  names  will 
be  explained  in  Appendix  B) . 

For  the  neighbor  regions  these  output  fluxes  become  inputs  and  have  to 
be  adjusted  according  to  the  differences  in  the  angle  DF1  of  the  different 
regions . 

The  simulation  starts  with  regions  in  the  sector  near  the 
breakdown  boundary.   At  this  time  there  is  no  data  for  input  flows  through 
faces  E  and  N  of  the  cell.   An  additional  run  of  the  whole  program  is  required 
in  order  to  take  these  calculated  flows  into  account.   If  the  accuracy  of  the 
results  is  important  we  may  run  this  type  of  iteration  several  times  until 
the  results  become  stable.   (Only  after  running  the  program  for  the  whole 
domain  once  we  shall  be  able  to  evaluate  the  importance  of  these  iterations.) 
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Figure  15.   Definition  of  input  and  output  flows  of  species  1  to  a  region  (KR) 
in  a  sector  (KS). 

(For  species  2  the  flux  names  will  change  as  follows: 
instead  FWN1(  )  *  FWN2(  ) 
instead  F0W1(  )  ♦  F0W2(  ) 
etc. 
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E.  COLLISIONLESS   FLOW 

In  several   sectors    near   the    breakdown   boundary  we   may   find   a  high   number 
density   and  the   mean  free   path   small   compared  with  the    size    of    a  cell.      There 
the   calculated  collisions    are   expected  to  have    an  influence    on  the    flow 
parameters.      For   wider   expansion   angles   the   collisions    become   rare   mainly 
because    of    the    decrease   in  the    density.      In  the    ambient    gas    the    mean  free    path 
(for    200   km  altitude)    is    240   m.      Comparing  this    number   with   the    size    of    the 
simulated   domain  may   lead  to  the   conclusion  that    there   the    flow  may   be 
regarded   as    collisionless. 

We    may   define    a  limiting   line    in  the   flow  where   the    collisions    become 
insignificant.      Consequently,    molecules    crossing  this    limit    will   in  fact 
continue    moving   in  straight    lines;    a  part    of    them  reach  the    solid  wall. 

Introducing  this    idea   of   the   collisionless    flow  we   may   reduce    the 
computation  time    and  the    memory   requirements. 

F.  TWO  DIMENSIONAL  PLANAR  FLOW  VS.  AXISYMMETRIC  FLOW 

The   cell   dimensions    are   completely   specified   only   in  three    dimensional 
flow.      When  applied  to  two   dimensional   flow,    the   simulation  may   be   regarded   as 
applying  to   an  arbitrarily   thin  slice    of    the   real   flow.      The   thickness    of   the 
slice    may   be   chosen   such   that    the    number    of   simulated   molecules    complies   with 
the   cell   volume    and   the    physical   number    density.      For    the    axisymmetric    flow  we 
have    defined  the    angle    DFI   as    the   third   coordinate   so  that    the    volume    of    the 
cell   is    completely   specified. 

Once   the    geometry   is    defined,    the    simulation   may   be    accomplished   and 
there    is    no   difference   if    doing   it    for    two   dimensional   or    for    axisymmetric 
flows. 
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APPENDIX  A  -  THE  AXSYM  PROGRAM 


A.l   DIFFERENT  REGIONS  IN  THE  JET 

For  the  two  dimensional  jet  with  initial  Mach  number  greater  than  unity 
the  different  regions  are  shown  in  Figure  (16), 
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Figure    16.     The   three   regions   in  an  underexpanded  jet. 

For   planar    2-D  flow  region   1   is   a  uniform  flow  core,    region  2  is    a  simple 
region  in  which  only  one   family  of   characteristics   define   the   flow,    and  region 
3  which  contains   the   intersection  of   the   two  families    of   characteristics. 
Because    our    intention  is    to   find   solutions    for    highly   underexpanded   jets   with 
very  low  ambient   pressure,    the   calculation  of   further   downstream  flow  is    not 
necessary. 

For   the   axisymmetric  ring  jet,   we   use   the   same   definition  for   the 
different   regions   however,   in  this   case   none   of   the   three   regions    has   uniform 
flow  and  is   not    a  simple   region. 
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As  shown  in  Equations  (9,10)  the  PM  function  (  v)  and  the  flow  direction 
(9)  of  a  point  at  location  I,j  may  be  calculated  from  the  v  and  9  of  two 
upstream  points  (I-l,j)  and  (I,j-1).   Later,  from  the  PM  function  at  the  new 
point  we  may  derive  the  local  Mach  number,  the  local  pressure,  temperature, 
velocity  and  other  thermodynamic  parameters  as  required. 

Definition  of  the  mesh  of  points  for  the  different  regions  is  shown  in 
Figure  (17). 

A. 2   PROGRAM  FLOWCHART 

A  simplified  flowchart  for  the  MOC  program  is  shown  in  Figure  (18).   The 
program  is  designed  to  solve  both  axisymmetric  as  well  as  two  dimensional  flow 

for  kD  =  2  it  solves  two  dimensional  flow 

for  kD  =  3  it  solves  axisymmetric  flow  (This  is  also  the  default 

condition) 
Initial  data  such  as  Mach  number  and  pressure  at  the  exit  surface,  ambient 
pressure  and  jet  gas  parameters  are  input  data. 

Output  data  contains  the  following  for  each  mesh  point: 

Mach  number,  coordinates  of  mesh  point  (R,X),  flow  direction  (TETA) , 

pressure,  temperature,  local  velocity,  Knudsen  number  basea  on  the 

distance  between  two  points  along  a  streamline,  mean  free  path  and 

breakdown  parameter  as  defined  by  Bird.* 

For  each  of  the  three  regions,  we  start  with  precalculated  boundary 
conditions  enabling  the  calculation  of  the  Mach  angles,  coordinates  of  mesh 
points  and  distances  d£  and  dn  as  described  in  [2J. 


*For  the  exit  plane  instead  the  Bird's  breakdown  parameter,  we  calculate  the 
ratio  between  time  per  three  collisions  and  time  of  motion.   Sometimes  this 
ratio  may  be  regarded  as  a  measure  of  the  breakdown  of  the  continuum  theory. 

40 


Figure    17.      Indexing   of   mesh  points   for    Che    different    regions    in  the    'AXSYM* 
program. 
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Figure    18.      AXSYM  Program  Flowchart 
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The  number  of  characteristics  used  in  the  program  is  arbitrary  and 
depends  on  the  required  resolution  (it  may  affect  also  the  accuracy  of  the 
results  and  the  amount  of  computation).   We  start  with  20  characteristics 
along  the  exit  cross  section  and  with  50  characteristics  in  the  Prandtl-Meyer 
fan  thus  a  total  of  70  characteristics  of  each  family  are  calculated.   In 
region  1  there  are  20  left  running  and  20  right  running  characteristics.   In 
region  2  there  are  20  right  running  and  50  left  running  characteristics.   In 
region  3  there  are  50  left  running  and  50  right  running  characteristics. 

In  region  3  we  limit  the  calculation  where  the  two  characteristics 
defining  a  new  mesh  point  intersect  at  an  angle  smaller  than  the 
computational  accuracy.   In  fact  this  occurs  far  downstream  where  continuum 
theory  becomes  invalid. 

After  defining  the  mesh  geometry  (successively)  we  calculate  the 
Prandtl-Meyer  function  and  flow  direction,  using  equations  (9,10). 

The  local  Mach  number  is  an  implicit  function  of  the  Pradtl-Meyer  angle. 
It  is  calculated  by  iterations  with  an  initial  guess  of  local  Mach  number 
set  equal  to  a  precalculated  Mach  number  at  an  adjacent  point,  and  the  slope 
of  the  function  as  given  by  Equation  (5).   Figure  (19)  shows  the  iterative 
procedure  for  evaluating  the  Mach  number  at  each  mesh  point. 
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Figure    19.      Iterative   procedure    for    Mach   number   calculation. 

Once    the    local   Mach   number   has    been  found   all   other    flow 
parameters    may   be    defined   using  the   Equation    (1,2,3). 
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Figure  20.  The  mash  of  characteristics  in  Region  2  and  Region  3 
Axisymmetric  ring  jet.   MQ  =  4.   Alt it ude= 200km. 
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Figure  21.   The  mesh  characteristics  in  Region  2  and  Region  3 
Axisymmetric  ring  jet.   MQ  =  2.   Altitude  ■  200km. 
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A. 3   PROGRAM  'AXSYII'  LISTING 

$J0B  AXS00010 

C   PROGRAM  AXSYM  AXS00020 

C  AXS00030 

!C  THIS  PROGRAM  CALCULATES  THE  ISENTROPIC  EXPANSION  OF  A  JET  BY  MEANS  OFAXS00040 

C   THE  METHOD  OF  CHARACTERISTICS.  AXS00050 

jC        FOR  A  TWO  DIMENSIONAL     JET  'KD'  SHOULD  BE  SET  EQUAL  TO  2  AXS00060 

iC        FOR  AN  AXISYMMETRIC  RING  JET  »KD»  SHOULD  BE  SET  EQUAL  TO  3  AXS00070 

jC  AXS00080 

IMPLICIT  REAL*8(A-H,0-Z,$)  AXS00090 

DIMENSION  TETA(20,50),AM(20,50),R(20,50),X(20,50),PM(20,50)  AXS0  0100 

DIMENSION  AMCOR(20,20),TETAC(20,20),XC(20,20),RC(20,20),PMC(20,20)AXS00110 

DIMENSION  DENSFC20,50)  AXS00120 

DIMENSION  AMX(50,50),TETAX(50,50),XX(50,50),RX(50,50),PMX(50,50)   AXS00130 

C  AXS0Q140 

C  TETA  IS  THE  FLOW  ANGLE  (RADIANS)  MEASURED  FROM  X  AXIS.  AXS00150 

C  AM  IS  THE  MACH  NUMBER  AXS00160 

C  R  IS  THE  RADIUS           (NORMAL  TO  THE  WALL)  AXS00170 

C  X  IS  THE  AXIAL  LOCATION    (PARALLEL  TO  THE  WALL)  AXS00180 

C  PM  IS  THE  PRANDTL  MEYER  FUNCTION  AXS00190 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0020  0 

C  THE  FOLLOWING  IS  DATA  FOR  THE  SPECIFIC  PROBLEM  AXS00210 

PAMB  =  8.4736E-5  AXS00220 

KD  =  3  AXS00230 

C    FOR   TWO    DIMENSIONAL    FLOW    KD=2    ^FOR    AXISYMMETRICAL    FLOW    KD=3  AXS00240 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0  0  250 

C  AXS00260 

C  CONSTANTS  AXS00270 

PI     =  3.141593  AXS00280 

BOLTZ  =  1.38032E-23  AXS00290 

AVOG   =  6.0225E+26  AXS00300 

RG     =  8314.3  AXS00310 

C  EXIT  SURFACE  AXS00320 

AMO  =  4.00  AXS00330 

TO   =  300.0  AXS00340 

PO   =  136.0  AXS00350 

C  AXS00360 

Rl   =  2.5  AXS00370 

XI   =  0.5  AXS00380 

C  AXS00390 

C  GAS  DATA  AXS00400 

GAMA  =  1.535  AXS00410 

DIAM  =  2.95E-10  AXS00420 

GM    =  17.0  AXS00430 

RJ  =  RG/GM  AXS00440 

CXS  =  PI*DIAM*DIAM  AXS00450 

GMM  =  GM/AVOG  AXS00460 

C  AXS00470 

C  MESH  DEFINITION  AXS00480 

C                 N1=CHARACTERISTICS  FROM  THE  EXIT  PLANE  AXS00490 

C                 N2=CHARACTERISTICS  FROM  THE  CORNERS  AXS00500 

Nl  =  20  AXS00510 

N2  =  50  AXS00520 

C  AXS00530 

C  CONSTANTS  FOR  THE  ISENTROPIC  EXPANSION  AXS00540 

Al  =  (GAMA-1.0VGAMA  AXS00550 

Bl  =  1.0/A1  AXS00560 

C  AXS00570 

A2  =DSQRT((GAMA+1.)/(GAMA-1. ))  AXS00580 

B2  =  1.0/A2  AXS00590 

C  AXS00600 

A3  =  (GAMA-D/2.  AXS00610 

C  AXS00620 

C     C  =  MACHXMACHXA3  +  1.  AXS00630 

C      D  =  MACHXMACH  -1 .  AXS00640 

CXXXXXX»X*XXXXXX»XXXXXX»XX*XXXXXXXXXXXXXXXXXXX»XXXXXKXXXXXXXXXXXXXXXX 

C  AXS00670 

C  DEFINE  STAGNATION  PARAMETERS  AXS00680 

C  AXS00690 

CXXXXXXXXXXKXXKXXXKKXXXXXXXXXXKXXXXXXXXXKXXXXXXXKXXXXXXXXKXXXXXXKXXXXXX  AXS00700 

C  =  (1.0+A3XAM0XAM0)  AXS00710 

C  PRESSURE  AXS00720 

ZL7 


PN  =  CXXB1  AXS0073C 

PSTG  =  POXPN  AXS0074C 

C  TEMPERATURE  AXS0075C 

TSTG  =  TOXC  AXS0076C 

C  DENSITY  AXS0077C 

DO  =  PO/CRJXTO)  AXS0078C 

DSTG  =  PSTG/CRJXTSTG)  AXS0079C 

C  AXS0080C 

Cxxxxxxxxxxxxx  AXS0081C 

WRITEC 6,1) PSTG, TSTG, DSTG, DIAM,GM  AXS0  082C 

1  FORMATCO', 'STAGNATION    PRESSURE= ■ , E12 . 5, ■  TEMPERATURE= ■ , FIO . 5, AXS0083C 
1'    DENSITY=',E12.5, »    MOL . DIAM= • , E12 . 5, ■  MOL .MASS= ' , FIO . 5)     AXS0084C 

WRITE  (6,2)P0,T0,D0,AM0  AXS0085C 

2  FORMAT('0','EXIT  PLANE    PRESSURE^ , E12 .5, «  TEMPERATURE= ' , FIO .5, AXS0086C 
1'    DENSITY=',E12.5, •    MACH= • , FIO . 5)  AXS0087C 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0  088C 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0089C 

C  AXS8090C 

C  DEFINE  FREE  STREAM  PARAMETERS  AT  THE  RIGHT  CORNER  AXS0091C 

C  AXS0092C 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0093C 

C  MACH  NUMBER  AXS0094C 

FSM  =  DSQRT(((PSTG/PAMB)XXA1-1.)/A3)  AXS0095C 

C  TEMPERATURE  AXS0096C 

FST  =  TSTG/C1.+A3XFSMXFSM)  AXS0097C 

C  DENSITY  AXS0098C 

FSD  =  PAMB/CRJXFST)  AXS0099C 

C  MACH  ANGLES  FOR  HEAD  AND  TAIL  OF  FAN  AXSOIOOC 

AMIT  =DARSIN(1./AM0)  AXS0101C 

AMIH  =DARSIN(1./FSM)  AXSO102C 

C  PRANDTL  MEYER  FUNCTION  FOR  HEAD  AND  TAIL  OF  FAN  AXS0103C 

Dl  =DSQRT(AM0XAM0-1.)  AXS0104C 

D2  =DSQRT(FSMXFSM-1.)  AXS0105C 

PMH  =  A2XDATAN(B2XD2)-DATAN(D2)  AXS0106C 

PMT  =  A2XDATAN(B2XD1)-DATAN(D1)  AXS0107C 

C  EXTERNAL  TURNING  ANGLE  (FREE  STREAM  ANGLE)=EXTA  AXS0108C 

EXTA  =  PMH  -  PMT  AXS0109C 

C  PRANDTL  MEYER  FAN  ANGLE  PMFA  AXS0110C 

PMFA  =  EXTA  -  AMIH  +  AMIT  AXS0111C 

C  AXS0112C 

Cxxxxxxxxxxxxxx  AXS0113C 

C      WRITE(6,3)PAMB,FST,FSD,FSM  AXS0114C 

C    3  FORMAT( '0', 'FREE  STREAM   PRESSURE= » , E12 . 5, »  TEMPERATURE= * , FIO . 5, AXS0115C 

C     1'    DENSITY=',E12.5, •    MACH= • , Fl 0 . 5///)  AXS0116C 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0117C 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0118C 

C  AXS0119C 

C  DEFINE  THE  MESH  OF  CHARACTERISTICS-AND  FLOW  PARAMETERS  AXS0120C 

C  AXS0121C 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0I22C 

C  THE  CORNER  POINT  -  LEFT  RUNNING  CHAR.  AXS0123C 

C  PMFA  IS  DEVIDED  INTO  N2  (=50)  NONEQUAL  DIVISIONS  AXS0124C 

C  AXS0125C 

RAT    =    (FSM/AMO-D/0.02  AXS0126C 

RAT2    =    FL0AKN2-1)  AXS0127C 

El  =  DL0G(RAT)/DL0G(RAT2)  AXS0128C 

C  AXSO1290 

WRITE(6,4)  AXSO13O0 

4  FORMAK »1', 'PRANDTL  MEYER  FAN           LINE  MACH         AXS01310 

1       P.M.  ANGLE         TETA'/)  AXS01320 

DO  20  N=2,N2  AXS01330 

EN2  =  FLOAKN-2)  AXS01340 

EN1  =  FLOAT(N-l)  AXS0135C 

AMB  =  AM0X(1.+.02X(EN2)XXE1)  AXS01360 

AMF  =  AMOX(1.+.02X(EN1)XXE1)  AXS01370 

DELM  =  AMF-AMB  AXS01380 

AM(1,N)  =  AMF  AXS01390 

Dl  =DSQRT(AM(1,N)XAM(1,N)-1. )  AXS01400 

PM(1,N)  =  A2XDATAN(B2XD1)-DATAN(D1)  AXS01410 

AMI  =DARSIN(1./AM(1,N))  AXS01420 

TETA(1,N)  =  PI/2.-(PM(l,N)-PMT)  AXS01430 

C  TETA  IS  THE  FLOW  ANGLE  ON  THE  CHARACTERISTICS  AT  THE  CORNER            AXS01440 


WRITE(6,10) 

10  FORMATC  0','    I 

J       MACH 

R 

X 

TETA 

1TEMP       PRESS 

VELOCITY 

KNUDSEN 

MFP 

ND 

1    P») 

c 

c 

R(1,N)  =  Rl  AXS01450 

X(1,N)  =  XI  AXS01460 

ALFAL  =  TETA(1,N)+AMI  AXS01470 

C  ALFAL  IS  THE  ANGLE  OF  THE  LEFT  RUNNING  CHARACTERISTICS  AXS01480 

C  AXS01490 

C  PRESSURE,  TEMPERATURE  AND  DENSITY  VARIATION  AT  THE  CORNER             AXS01500 

C  =  AM(1,N)*AM(1,N)XA3+1.  AXS01510 

PRES  =  PSTG/(CXXB1)  AXS01520 

TEMP  =  TSTG/C  AXS01530 

DENSF(1,N)  =  PRES/CRJXTEMP)  AXS01540 

C  AXS01550 

Cxxxxxxxx  AXS01560 

WRITE  (6,5)  N,AM(1,N),PM(1,N),TETA(1,N)  AXS01570 

5  FORMATC  ',»  ' , 115, 3E20 .5)                AXS01580 

Cxxxxxxxx  AXS01590 

20  CONTINUE  AXS01600 

C  AXS01610 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS01620 

C  AXSD1630 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS01640 

AXS01650 
AXS01660 
AXS01670 
AXS01680 
AXS01690 
AXS01700 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS01710 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS017  20 
C  AXS01730 

C  DEFINE  THE  PARAMETERS  AT  THE  EXIT  SURFACE  (PLANE)  AXS01740 

C  AXS01750 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS 01760 

DO  25  J  =  1,N1  AXS01770 

AMCORCJ)  =  AMO  AXS01780 

D  =  DSQRT(AMOXAMO-L)  AXS01790 

PMC(1,J)  =  A2XDATAN(B2XD)-DATAN(D)  AXS01800 

25  TETAC(1,J)  =  PI/2.  AXS01810 

C  AXS01820 

C  THE  EXIT  PLANE  IS  DIVIDED  INTO  (MM)  DIVISIONS  (Nl  POINTS)  AXS01830 

DO  30  J=l,20  AXS01840 

SOUND  =  DSQRT(GAMAXRJXTO)  AXS01850 

VELO  =  AMCORd,  J)XSOUND  AXS01860 

DNO  =  P0/(BOLTZXT0)  AXS01870 

FPO  =  .707/(DN0xCXS)  AXS01880 

DX  =  X1X2./FL0AT(N1)  AXS01890 

AKNO  =  FPO/DX  AXS01900 

C  AXS01910 

XC(1,J)  =  X1X(1.-FL0AT(J-1)/FL0AT(N1-1)X2.  )  AXS01920 

CENTR  =  DABS(XC(1,J))  AXS01930 

IF(CENTR.LT. 0.001)  XC(1,J)  =  0.  AXS019^0 

RC(1,J)  =  Rl  AXS01950 

C  AXS01960 

C AXS01970 

C     KZ=lxxxxxxxxxxxxxxxxxxx  AXS01980 

IF  (J.GT.l)  GO  TO  29  AXS01990 

C  AT  THE  EXIT  PLANE  THE  BREAKDOWN  PARAMETER  IS  EVALUATED  BY  MEANS  OF     AXS02000 

C  THE  RATIO  OF  THE  COLLISION  TIME  AND  THE  FLOW  TIME.  LENGTH  SCALE  IS     AXS02010 

C  THE  MESH  DIMENSION.  AX502020 

SCALE  =  DXxDTAN(DARSIN(l./AM0))/2.  AXS02030 

TIME1  =  SCALE/VELO  AXS02040 

TIME2  =  3./(4.xCXSXDN0xDSQRT(BOLTZXT0/(PIXGMM)))  AXS02050 

P  =TIME2/TIME1  AXS02060 

DENSF(1,1)  =  DO  AXS02070 

29  CONTINUE  AXS02080 
C  AXS02090 

WRITE(6,ll)l,J,AMCOR(l,J),Rl,XC(l,J),TETAC(l,J),T0,P0,VEL0,AKN0,FPAXS02100 

10, DNO, P  AXS02110 

11  FORMATC  »,2I4,5F10.3,3E12.3,F9.4,2E12.3)  AXS02120 

30  CONTINUE  AXS02130 
C  AXS02140 

CXXX*XXXXXXXXXXXXX*XXXXX*XXXXXXXXXXXXXXXXXXXXXXX**XXXXXXXXXX*XXXXXXXXXXXAXS02150 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0216  0 

'♦9 


C  AXS0217 

C  CALCULATE  THE  FLOW  PARAMETERS  IN  THE  CORE  BOUNDED  BY  THE  TWO  MACH  AXS0218 

C  WAVES  STARTING  AT  THE  NOZZLE  LIPS  (CORNER  POINTS)  AXS0219) 

C  AT  THE  EXIT  THE  FLOW  IS  ASSUMED  TO  BE  UNIFORM  AXS0220 

C  AXS0221 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0222 

WRITE  (6,198)  AXS0223 

198  FORMATCl','    CORE    »/)  AXS0224 

APM=0.  AXS0225 

BPM=0.  AXS0226 

CPM=0.  AXS0227 

C  AXS0228 

DO  199  I  =  2,N1  AXS0229 

WRITE  (6,10)  AXS0230 

DO  199  J  =  1,N1  AXS0231 

IF  ((I  +  J-D.GT.N1)  GO  TO  199  AXS0232 

AMIL  =  DARSIN(l./AMCOR(I-l,J))  AXS0233 

ALFAL  =  PI-(TETAC(I-1,J)+AMIL)  AXS0234 

C  AXS0235 

AMIR  =  DARSIN(l./AMCOR(I-l,J  +  D)  AXS0236 

ALFAR  =  TETAC(I-1,J+1)-AMIR  AXS0237 

C  AXS0238 

XC(I,J)  =  (RC(I-1, J)-RC(I-1,J+1)+XC(I-1,J)XDTAN(ALFAL)+XC(I-1, J+DAXS0239 

1*DTAN(AL FAR) )/( DTAN( ALFAL )+DTAN( ALFAR))  AXS0240 

CENTR  =  DABS(XC(I,J))  AXS0241 

IF  (CENTR. LT. 0.001)  XC(I,J)  =  0.  AXS0242 

RC(I,J)  =  RC(I-1, J+1)+(XC(I, J)-XC(I-1,J+1))*DTAN(ALFAR)  AXS0243 

C  AXS0244 

DKSI  =  DSQRT((XC(I,J)-XC(I-1,J  +  1))**2+(RC(I,J)-RC(I-1,J  +  1))X*2)    AXS0245 

DETA  =  DSQRT((XC(I,J)-XC(I-1,J))XX2+(RC(I,J)-RC(I-1,J))*X2)  AXS0246 

C  AXS0247 

C  AXS0248 

C  CALCULATE  NOW  THE  PRANDTL  MEYER  FUNCTION  AND  FLOW  ANGLE  IN  CORE.  AXS0249 

APM  =  PMC(I-1,J)+PMC(I-1,J+1)+TETAC(I-1,J+1)-TETAC(I-1,J)  AXS0250 

IF  (KD.EQ.2)  GO  TO  151  AXS0251 

BPM  =  DSIN(AMIR)XDSIN(TETAC(I-1,J+1))/RC(I-1,J+1)XDKSI  AXS0252 

CPM  =  DSIN(AMI'L)XDSIN(TETAC(I-1,J))/RC(I-1,J)XDETA  AXS0253 

151  PMC(I,J)  =  (APM+BPM+CPM)/2.0  AXS0254 

C  AXS0255 

APM  =  PMC(I-1,J+1)-PMC(I-1,J)+TETAC(I-1,J+1)+TETAC(I-1,J)  AXS0256 

TETAC(I,J)  =  (APM+BPM-CPM)/2.0  AXS0257 

C  AXS0258 

C  DEFINE  MACH  NUMBER  (BY  ITERATIONS)  AXS0259 

C  INITIAL  GUESS  AMCOR( I , J ) =AMCOR( 1-1 , J+l )  AXS0260 

C  AXS0261 

AMG  =  AMCOR(I-l,J+l)  AXS0262 

KZ  =  0  AXS0263 

154  IF  (KZ.GE.100)  GO  TO  160  AXS0264 

KZ  =  KZ+1  AXS0265 

C  =  AMGXAMGXA3+1.  AXS0266 

D  =  DSQRT(AMGXAMG-1.)  AXS0267 

PMCAL  =  A2XDATAN(B2XD)-DATAN(D)  AXS0268 

DELNI  =  PMCAL  -  PMC(I,J)  AXS0269 

DEL  =  DABS(DELNI)  AXS0270 

IF  (DEL. LT. .000002)  GO  TO  160  AXS0271 

IF  (DELNI. LT. 0. )  GO  TO  156  AXS0272 

AMG  =  AMGX.999  AXS0273i 

GO  TO  154  AXS0274I 

156  AMG  =  AMGX(1.-DELNIXC/D)  AXS0275I 

GO  TO  154  AXS0276I 

160  AMCOR(I,J)  =  AMG  AXS0277I 

C AXS0  27  8I 

C  AXS0279I 

C  CALCULATE  FLOW  PARAMETERS  AXS0280I 

IF  (J.GT.l)  GO  TO  197  AXS0281I 

C  =  AMC0R(I,J)XAMC0R(I,J)XA3+1.  AXS0282I 

PRES  =  PSTG/(CXXB1)  AXS0283I 

TEMP  =  TSTG/C  AXS0284I 

DN    =  PRES/(B0LTZXTEMP)  AXS0285I 

FP    =  .707/(DNXCXS)  AXS0286I 

SCALE  =  DKSI  x  DSIN(ALFAL)  AXS0287I 

AKN  =  FP/SCALE  AXS0288I 
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SOUND  =  DSQRTCGAMAXRJXTEMP)  AXS02890 

VEL  =  SOUND*AMCOR(I,J)  AXS02900 

C  AXS02910 

C  BREAKDOWN  PARAMETER  AS  DEFINED  BY  'BIRD*.  AXS02920 

DENSF(I,J)  =  PRES/CRJXTEMP)  AXS02930 

DDENS  =  DENSF(I-1,J)  -  DENSF(I,J)  AXS02940 

COLF   =  4.XCXSXDN*DSQRT(B0LTZXTEMP/(PI*GMM))  AXS02950 

P  =  VELXDDENS/(SCALEXDENSF(I,J)xCOLF)  AXS02960 

197  CONTINUE  AXS02970 

C  AXS02980 

C      TIME1  =  SCALE/VEL  AXS02990 

C     TIME2  =  3./(4.XCXSXDNXDSQRT(B0LTZXTEMP/(PlxGMM)))  AXS03000 

C  P  =  TIME2/TIME1  AXS03010 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS 03020 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS03030 

Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*AXS0  304  0 

WRITE(6,11)I,J,AMCOR(I,J),RC(I,J),XC(I,J),TETAC(I,J),TEMP,PRES,VELAXS03050 

1,AKN,FP,DN,P  AXS03060 

199  CONTINUE  AXSJ33070 
C  AXS03080 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0  309  0 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0  310  0 

C  MATCH  CORE  AND  FAN  POINTS  AXS03110 

DO  200  1=1, Nl  AXS03120 

X(I,1)  =XC(I,1)  AXS03130 

R(I,1)=  RC(I,1)  AXS03140 

AM(I,1)  =  AMCOR(I,l)  AXS03150 

PM(I,1)  =  PMC(I,1)  AXS03160 

200  TETA(I,1)  =  TETAC(I,1)  AXS03170 
C  AXS03180 

CXXXXXXXX*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX*XX*AXS0319  0 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0320  0 

C  AX503210 

C  CALCULATE  FLOW  PARAMETERS  IN  REGION  2  (SIMPLE  PRANDTL  MEYER  FAN).      AXS03220 

C  AXS03230 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0  3240 

NRITE  (6,298)  AXS03250 

298  FORMAT  (U','      REGION     2*/)  AXS03260 

C  AXS03270 

KFIN  =  51  AXS03280 

C  AXS03290 

DO  299  I  =  2,N1  AXS03300 

NRITE(6,10)  AXS03310 

DO  299  J  =  2,N2  AXS03320 

IF  (J.GT.KFIN)  GO  TO  299  AXS03330 

KZ  =  0  AXS03340 

AMIL  =DARSIN(1./AM(I-1, J))  AXS03350 

AMIR  =DARSIN(1./AM(I,  J-D)  AXS03360 

ALFAL  =  PI-(TETA(I-1,J)+AMIL)  AXS03370 

ALFAR  =  TETA(I,J-1)-AMIR  AXS03380 

C  AXS03390 

C  CHECK  ANGLES  AND  CALCULATE  COORDINATES  AXS03400 

ANGLE1  =  PI/2.-. 000001  AXS03410 

ANGLE2  =  PI/2. +.000001  AX503420 

IF  (ALFAL. LE.ANGLE1. OR. ALFAL. GE.ANGLE2)  GO  TO  201  AX303430 

X(I,J)  =  X(I-1,J)  AXS03440 

GO  TO  207  AXS03450 

201  IF  (ALFAR. LE.ANGLE1. OR. ALFAR. GE.ANGLE2)  GO  TO  205  AXS03460 
X(I,J)  =  X(I,J-1)  AXS03470 
GO  TO  207  AXS03480 

C  AXS03490 

205  X(I,J)  =  (R(I-1,J)-R(I,J-1)+X(I,J-1)XDTAN(ALFAR)+X(I-1, J )XDTAN( ALFAXS03500 

1AL))/(DTAN(ALFAL)+DTAN( ALFAR))  AXS0  3510 

C  AXS03520 

207  IF  (ALFAL. LE. 0.000001. OR. ALFAL. GE. 0.000001)  GO  TO  209             AXS03530 

R(I,J)  =  R(I-1,J)  AXS03540 

GO  TO  213  AXS03550 

209  IF  (ALFAR. LE. 0.000001. OR. ALFAR. GE. 0.000001)  GO  TO  211             AXS03560 

R(I,J)  =  R(I,J-1)  AX503570 

GO  TO  213  AXS03580 

C  AXS03590 

211  R(I,J)  =  (X(I, J)-X(I, J-1))XDTAN(ALFAR)+R(I,J-1)  AXS03600 
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213  DKSI  =  DSQRT((R(I,J)-R(I,J-1))XX2+(X(I,J)-X(I,J-1))XX2) 
DETA  =  DSQRT(CRCI,J)-R(I-1,J))*X2+(X(I,J)-XCI-1,J))XX2) 
C 

IF  CRCI,J).GT.O..AND.X(I,J).GT.O.)  GO  TO  219 

KFIN  =  J-l 

WRITE(6,12) 
12  FORMATC  ' ,»  FURTHER  POINTS  ON  THE  CHARACTERISTICS  ARE  DIVERGENT' 
219  CONTINUE 
C  LOCATION  OF  THE  NEW  MESH  POINT  HAS  BEEN  FOUND 
C 
C 
C  CALCULATE  NOW  PRANDTL  MEYER  FUNCTION  AND  FLOW  DIRECTION  FOR  NEW  POINT 

APM  =  PM(I,J-1)+PM(I-1,J)-TETA(I-1,J)+TETA(I,J-1) 

IF  (KD.EQ.2)  GO  TO  251 

BPM  =  DSIN(AMIR)*DSIN(TETA(I,J-1))/R(I,J-1)XDKSI 

CPM  =  DSIN(AMIL)XDSIN(TETA(I-1,J))/R(I-1,J)XDETA 
251  PM(I,J)  =  .5x(APM+BPM+CPM) 
C 

APM  =  PM(I,J-1)-PM(I-1,J)+TETA(I,J-1)+TETA(I-1,J) 

TETA(I,J)  =  .5X(APM+BPM-CPM) 
C 

C  CALCULATE  NOW  MACH  NUMBER  FOR  EACH  POINT 
C  INITIAL  GUESS  AM(I,J)  =  AMCI-1,J) 

AMG  =  AM(I-1,J) 

KZ  =  0 
254  IF  (AMG. GT. 200.)  GO  TO  257 

D  =  DSQRTCAMGXAMG-1.) 

GO  TO  258 

257  D  =  AMG 

258  IF  (KZ.GE.100)  GO  TO  260 
KZ  =  KZ  +  1 

PMCAL  =  A2XDATAN(B2XD)-DATAN(D) 
•   DELNI  =  PMCAL  -  PM(I,J) 
DEL  =  DABS(DELNI) 
IF  (DEL. LT. .000002)  GO  TO  260 
IF  (DELNI. LT.O. )GO  TO  256 
AMG  =  AMGX.999 
GO  TO  254 

c////////////////////////////////////////////////////////////////////// 

256  IF  (AMG. LT. 2000. )GO  TO  2560 
KFIN  =  J-l 
GO  TO  299 

c////////////////////////////////////////////////////////////////////// 

2560  Dl  =  (A3XAMGXAMG+1 . ) 

AMG    =    AMGXU.-DELNIXD1/D) 

GO    TO    254 
C 

260  AM(I,J)  =  AMG 
C  CALCULATE  NOW  LOCAL  TEMPERATURE, PRESSURE, VELOCITY, KNUDSEN  NO, 

C  =  AM(I,J)XAM(I, J)XA3+1 

PRES  =  PSTG/(CXXB1) 

TEMP  =  TSTG/C 

DN  =  PRES/(BOLTZXTEMP) 

DENSF(IJ)  =  PRES/(RJXTEMP) 

DDENS  =  DENSF(I-1, J-l )-DENSF( I , J) 

SCALE  =  DSORT((X(I,J)-X(I-l, J-l )  )*x2  +  ( R( I , J )-R( 1-1 , J-l ) )XX2) 
Cxxxxxx*xxxxxx*x**xx*xxxx*x**xxx*xxx*xxxx**x 

c 

271  FP  =  .707/(DNxCXS) 
C 

AKN  =  FP/SCALE 

SOUND  =DSQRT(GAMAXRJXTEMP) 

VEL  =  AM(I, J)XS0UND 


C 

C 


COLF  =  4.XCXS*DNXDSQRT(B0LTZXTEMP/(PIXGMM)) 
P  =  VELXDDENS/(SCALEXDENSF(I,J)XC0LF) 


C  PRINT  RESULTS  FOR  MESH  POINTS 
KL  =  (-l)xxi 
IF  (KL.LT.O)  GO  TO  299 
WRITE(6,11)I, J,AM(I, J),R(I, J  ) , X( I , J  ) , TETA( I , J ) , TEMP , PRES, VEL , AKN, 


AXS036K 
AXS0362( 
AXS03631 
AXS0364( 
AXS0365( 
AXS0366( 

)AXS0367( 
AXS0368( 
AXS0369I 
AXS0370I 
AXS037K 
AXS0372I 
AXS0373I 
AXS0374I 
AXS0375I 
AXS0376I 
AXS0377I 
AXSJ3378I 
AXS0379I 
AXS0330I 
AXS0381I 
AXS0382I 
AXS0383I 
AXS0384I 
AXS0385I 
AXS0386I 
AXS0387I 
AXS0388I 
AXS0389I 
AXS0390I 
AXS0391I 
AXS0392I 
AXS0393I 
AXS0394I 
AXS0395I 
AXS0396I 
AXS0397I 
AXS0398I 
AXS0399I 
AXS0400I 
AXS0401I 
AXS0402I 

/AXS0403I 
AXS0404I 
AXS0405I 
AXS0406I 
AXS0407I 
AXS0408I 
AXS0409I 
AXS0410I 
AXS0411I 
AXS0412I 
AXS0413I 
AXS0414I 
AXS0415I 
AXS0416I 
AXS04171 
AXS04181 
AXS0419( 
AXS0420C 
AXS042K 
AXS0422C 
AXS0423C 
AXS0424C 
AXS0425C 
AXS0426C 
AXS0427C 
AXS0428C 
AXS0429C 
AXS0430C 
AXS0431C 
AXS0432C 


1FP,DN,P  AXS04330 

299  CONTINUE  AXS04340 
C  AXS04350 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0436  0 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS 04370 

C  AXS04380 

C  MATCH  'REGION  2'  AND  'REGION  3'  POINTS  AXS04390 

C  AXS04400 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS04410 

L  =  KFIN-  1  AXS04420 

DO  300  J  =  1,L  AXS04430 

XX(1,J)     =  X(20,J)  AXS04440 

RX(1,J)     =  R(20,J)  AXS04450 

AMX(1,J)    =  AM(20,J)  AXS04460 

PMX(1,J)    =  PM(20,J)  AXS04470 

300  TETAX(1,J)  =  TETA(20,J)  AXS04480 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS 0449  0 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0450  0 

C  AXS04510 

C  CACULATE  FLOW  PARAMETERS  FOR  REGION  3  AXS04520 

C  AXS04530 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS0454  0 

WRITE  (6,397)  AXS04550 

397  FORMAT  (•IS1    REGION     3V)  AXS04560 

DO  399  I  =  2,L  AXS04570 

DO  399  J  =  1,1  AXS04580 

IF  (J.GT.KFIN)  GO  TO  399  AXS04590 

IF  (J.GT.I)GO  TO  320  AXS04600 

WRITE  (6,10)  AXS04610 

320  KZ  =  0  AXS04620 

C  AXS04630 

AMIL  =  DARSIN(1./AMX(I-1,J))  AXS04640 

ALFAL  =  PI-(TETAX(I-1,J)+AMIL)  AXS04650 

IF  (J.GT.I)  GO  TO  301  AXS04660 

ALFAR  =  ALFAL  AXS04670 

TETAX(I,J)  =  PIX.5  AXS04680 

TETAX(I,J-1)  =  PI-TETAX(I-1,J)  AXS04690 

XX(I,J)  =  0.  AXS04700 

RX(I,J)  =  RX(I-1,J)+XX(I-1, J)XOTAN(ALFAL)  AXS04710 

RX(I,J-1)  =  RX(I-1,J)  AXS04720 

XX(I,J-1)  =  -XX(I-1,J)  AXS04730 

DKSI  =  (RX(I,J)-RX(I-1,J))/DSIN(ALFAL)  AXS04740 

•  DETA  =  DKSI  AXS04750 

PMX(I,J-1)  =  PMX(I-1,J)  AXS04760 

GO  TO  316  AXS04770 

301  AMIR  =  DARSIN(1./AMX(I,J-D)  AXS04780 
ALFAR  =  TETAX(I,J-1)-AMIR  AXS04790 

C  AXS04800 

IF  (ALFAL. LE.ANGLE1. OR. ALFAL. GE.ANGLE2)  GO  TO  302  AX504810 

XX(I,J)  =  XX(I-1,J)  AXS04820 

GO  TO  307  AXS04830 

302  IFCALFAR.LE.ANGLE1. OR. ALFAR. GE.ANGLE2)  GO  TO  305  AXS04840 
XX(I,J)  =  XX(I,J-1)  AXS04350 
GO  TO  307  AXS04860 

305  XX(I,J)  =  (RX(I-1, J)-RX(I, J-1)+XX(I, J-1)XDTAN(ALFAR)+XX(I-1,J)XDTAAXS04870 

1N( AL FAL ) )/(DTAN(AL FAD +DTAN( ALFAR))  AXS 0438  0 

C  AXS04890 

307  IF  (ALFAL. LE. 0.000001. OR. ALFAL. GE. 0.000001)  GO  TO  309             AXS04900 

RX(I,J)  =  RX(I-1,J)  AXS04910 

GO  TO  315  AXS04920 

309  IF  (ALFAR. LE. 0.000001. OR. ALFAR. GE. 0.000001)  GO  TO  311             AXS04930 

RX(I,J)  =  RX(I,J-1)  AXS04940 

GO  TO  315  AXS04950 

311  RX(I,J)  =  RX(I,J-1)+(XX(I,J)-XX(I,J-1))XDTAN(ALFAR)  AXS04960 

C  AXS04970 

315  DKSI  =  DSQRT((RX(I,J)-RX(I,J-1))X*2+(XX(I,J)-XX(I,J-1))XX2)  AXS04980 
DETA  =  DSQRT((RX(I, J)-RX(I-1,J))XX2+(XX(I, J)-XX(I-1,J))XX2)  AXS04990 

316  CONTINUE  AXS05000 
C  AXS05010 

IF  (RX(I,J) .GE.O. .AND.XX(I,J) .GE.O. )GO  TO  319  AXS05020 

KFIN  =  J-l  AXS05030 

WRITE  (6,398)  AXS05040 
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398  FORMAT  ('  ',* FURTHER  POINTS  ON  THE  CHARACTERISTICS  ARE  DIVERGENT' 
C 

C  LOCATION  OF  THE  NEW  POINT  HAS  BEEN  FOUND 
CXXXXXXXXXXXXXXXXXXXXXX*XXXXXXXXXXXXXXXXXXXXXXXXXXX*XXXXXXXXXXXXXXXXXX 

c 

C  CALCULATE  NON  P.M.  ANGLE  AND  FLOW  DIRECTION 

319  APM  =  PMX(I,J-1)+PMX(I-1,J)-TETAX(I-1,J)+TETAX(I,J-1) 
IF(KD.EQ.2)  GO  TO  351 
BPM  =  DSIN(AMIR)xDSIN(TETAX(I, J-l ) )/RX( I , J-l )XDKSI 

CPM  =  DSIN(AMIL)*DSIN(TETAXCI-1,J))/RXCI-1,J)*DETA 
351  PMX(I,J)  =  .5x(APM+BPM+CPM) 

APM  =  PMX(I,J-1)-PMX(I-1,J)+TETAX(I,J-1)+TETAXCI-1,J) 

TETAX(I,J)  =  .5*(APM+BPM+CPM) 
C 

C  CALCULATE  NOW  THE  MACH  NUMBER  FOR  EACH  POINT 
C  INITIAL  GUESS  AMX(I,J)  =  AMX(I-1,J) 

AMG  =  AMX(I-1,J) 

KZ  =  0 

KH  =  0 

KL  =  0 
354  IF  (AMG. GT. 2000.)  GO  TO  357 

D  =  DSQRTCAMGXAMG-1.) 

GO  TO  358 

357  D  =  AMG 

358  IF  (KZ.GE.50)G0  TO  360 
KZ  =  KZ  +1 

PMCAL  =  A2*DATANCB2XD)-DATAN(D) 
DELNI  =  PMCAL-PMX(I,J) 
DEL  =  DABS(DELNI) 
IF  (DEL. LT. .000002)  GO  TO  360 
IF  (DELNI. LT.0.)GO  TO  356 
AMG  =  AMGX.98 
GO  TO  354 
356  IF  (AMG. LT. 5000.)  GO  TO  3560 
KFIN  =  J-l 
GO  TO  399 
3560  Dl  =  A3XAMGXAMG+1. 
Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
DDELNI  =  DELNI 
IF  (AMG.GT.20.)  DDELNI  =  DELNIX( . 95XXKZ) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

AMG  =  AMGX(1.-DDELNI*D1/D) 

GO  TO  354 
360  AMX(I,J)  =  AMG 
C 
C  CALCULATE  NOW  THE  LOCAL  TEMPERATURE,  PRESSURE,  VELOCITY, KNUDSEN 

C  =  AMX(I,J)XAMX(I, J)XA3+1. 

PRES  =  PSTG/(CXXB1) 

TEMP  =  TSTG/C 

DN  =  PRES/(B0LTZXTEMP) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

DENSF(I,J)=PRES/(RJXTEMP) 

DDENS  =  DENSF(I-1, J-l )-DENSF( I ,  J ) 

FP  =  .707/(DNXCXS) 

SCALE=DSQRT((XX(I,J)-XX(I-1,J-1))XX2+(RX(I,J)-RX(I-1,J-1))XX2) 

AKN=FP/SCALE 
C 
C 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C  ALFALL  =  -ALFAL+PI 

C  AF  =  ALFALL-PI/2. 

C  AF  =  DABS(AF) 

C  IF  (AF.LT. .000001)  GO  TO  370 

C  BF  =  l./DSQRTU.  +  DTAN(  AL  FALL  )x*2) 

C  SCALE=BFX(RX(I,J-1)-DTAN(ALFALL)X(XX(I,J-1)-XX(I-1,J))-RX(I-1,J)) 

C  SCALE  =  DABS(SCALE) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

SOUND  =  DSQRT(GAMAXRJ*TEMP) 
VEL  =  AMX(I,J)XS0UND 
C 

COLF=4.XCXSXDNXDSORT(BOLTZXTEMP/(PIXGMM)) 
P=VELXDDENS/(SCALEXDENSF(I,J)xCOLF) 


r)4 


C  AXS05770 

C  PRINT  RESULTS  AXS05780 

WRITE(6,  11)I,J,AMX(I,J),RXCI,J),XX(I,J),TETAXCI,J),TEMP,PRES,VEL,AXS05790 

1AKN,FP,DN,P  AXS05800 

399  CONTINUE  AXS05810 

STOP  AXS05820 

END  AXS05830 

SENTRY  AXS05840 
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A. 4   LIST  OF  MAIN  SYMBOLS  USED  IN  'AXSYM* . 


Parameter 

Name 

Physical   Name 

Units 

Type 

Description 

PAMB 

ambient    pressure 

pascals 

real 

ambient    atmosphere 
pressure 

KD 

K  dimensions 

integer 

KD=2   for    two   dimen- 
sional  jet 

KD=3   for    axisymmetric 
ring  jet 

PI 

IT 

"™ 

real 
constant 

BOLTZ 

Boltzmann   constant 

Joule 
degree 

real 

1.38032      10-" 
joules /degree 

AVOG 

Avogadro's    constant 

molecules 

real 

6.0225   x    10Zb 
1/kmol 

mol 

RG 

Universal   gas 
constant 

Joule 

real 

8314.3 

mol.deg 

AM£ 

Mach   No.    (E.P) 

~ 

real 

Mach   number    at    exit 
plane 

TO 

Temperature    (E.P) 

°k 

real 

Temperature    at    exist 
plane 

PO 

Pressure    (E.P) 

pascals 

real 

Pressure    at    exit 
plane 

Rl 

Cylinder   radius 

m 

real 

Radius    of    the 
cylindrical   vehicle 

XI 

0.5*   nozzle   width 

m 

real 

Half   width    of    nozzle 

GAMA 

"™ 

real 

averaged   heat 
capacity   ratio    (jet) 

DIAM 

Molecule    diameter 

m 

real 

averaged   molecular 
diameter    (jet) 

GM 

Molecular   mass 

kg/kmol 

real 

averaged   molecular 
weight    (jet) 

RJ 

Gas    constant 

joules 
kj.deg 

real 

gas    constant    (jet) 

CXS 

m^ 

real 

collision   cross 
section    (hard   sphere) 

GMM 

mass    of    a  molecule 

kg 

real 

averaged   mass    of    a 
molecule 

Nl 

integer 

number    of    divisions 
(characteristics    from 
the   exist    plane) 

N2 

integer 

number    of    character- 
istics   in  the    Prandtl 

Meyer    fan 

Al 

real 

Y  "    1 
Y 

Bl 

real 

J 

Y 

T-   1 
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A. 4   LIST  OF  MAIN  SYMBOLS  USED  IN  »AXSYM»  (CONTINUED) 


Parameter 
Name 

Physical  Name 

Units 

Type 

Description 

A2 

real 

Ky  +  D/(y  -  D] 

B2 

real 

1/A2 

A3 

real 

T-  1 
2 

C 

real 

M2  y  ~  1  +  1 
2 

D,D1,D2 

real 

li2-  -    1 

PSTG 

pascals 

real 

stagnation  pressure 

TSTG 

°k 

real 

stagnation  tempera- 
ture 

DSTG 

kg/m3 

real 

stagnation  density 

FSM 

real 

free  stream  Mach 
number 

FST 

°k 

real 

free  stream  tempera- 
ture 

FSD 

kg/m3 

real 

free  stream  density 

AMIT 

MT 

radius 

real 

Mach  angle  (tail  of 
P. 11.  fan) 

AMIH 

m 

radius 

real 

ilach  angle  (head  of 
P.M.  fan) 

PMT 

v^ 

radius 

real 

P.M.  function  (tail) 

PMH 

VH 

radius 

real 

P.M.  function  (tiead) 

RAT , RAT2 , 
El 

real 

used  to  define  a 
loparithmic  division 
of  the  corner  char- 
acteristics (50  lines 
in  P.M.  fan)(a  linear 
division  v/ould  have 
resulted  in  concen- 
tration of  character- 
istics at  high  Mach 
numbers ) 

DELM 

real 

difference  of  Mach 
numbers 

AM(I,J) 

11 

real  2-D 
array 

Mach  number  at 
location  (I,j)(used 
for  the  corner) 

PM(I,J) 

V 

radius 

P.M.  function  at 
location  (I,j)(used 
for  the  corner) 

AMI 

M 

radius 

real 

Mach  angle 

R(I,J) 

m 

real  2-D 
array 

radius  (or  ordinate) 
at  (I,j) 

TETA(I,J)    9 

radius 

real  2-D 
array 

flow  direction  at 
(I,j) 
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A. 4      LIST   OF   1IAIN   SYMBOLS   USED    IN    'AXSYM'    (CONTINUED) 


Parameter 

Name 

Physical  Name 

Units 

Type 

Description 

ALFAL 

radians 

real 

angle  of  a  left  run- 
ning characteristics 

ALFAR 

radians 

real 

angle  of  a  right  run- 
ning characteristics 

DENSF(I,J) 

kg/m3 

real  2-D 
array 

gas  density  at  point 
IJ 

AMCOR(I,J) 

real  2-D 
array 

Mach  number  at  mesh 
points  in  region  (1) 
(region  1  is  bounded 
by  the  two  tail 
characteris  tics 

TETAC(I.J) 

radians 

real  2-D 
array 

flow  direction  at 
mesh  points  (region 
1) 

PMC(I,J) 

radians 

real  2-D 
array 

P.M.  function  at  mesh 
points  (region  1) 

XC(I,J) 

m 

real  2-D 
array 

X  coordinate  at  mesh 
points  (region  1) 

RC(I,J) 

m 

real  2-D 
array 

radius  (ordinate)  at 
mesh  points  (region 
1) 

AMIL 

radians 

real 

Mach  angle  far  lert 
running  character- 
istics 

AMIR 

radians 

real 

Mach  angle  far  right 
running  character- 
is  tics 

DKSI 

<U 

m 

real 

distance  between  mesh 
points 

DETA 

dn 

m 

real 

distance  between  mesh 
points 

DELNI 

dv 

radians 

real 

P.M.  differential 

AMG 

radians 

Mach  number  (used  for 
iterative  calcula- 
tion) 

PRES 

pressure 

TEMP 

temperature 

DN 

number  density 

FP 

mean  free  path 

AKN 

Knudsen  number 

SOUND 

speed  of  sound 

VEL 

absolute  local 
velocity 

DDENS 

density  difference 
between  two  points  a- 
long  a  steamline 

COLF 

collision  frequency 

P 

breakdown  parameter 
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A. 5      AXSYM   PROGRAM  USER'S   GUIDE 

1.  Input  data: 

ambient  pressure  PAMB 

Mach  number  (Exit  plane)  AMO 

Temperature  (Exit  plane)  TO 

Pressure  (Exit  plane)  PO 

Half  width  of  nozzle  XI 

Radius  of  nozzle  ring  Rl 

Specific  heat  ratio  of  jet  gas  GAIIA 

average  molecular  diameter  (jet)  DIAM 

average  molecular  weight  Gil 

2.  Options  for  flow  geometry 

two  dimensional  flow  KD=2 

ring  jet  KD=3 

default  condition  KD=3 

3.  Resolution  of  mesh  points 

to  change  the  resolution  of  the  mesh  points 
a  -  change  Nl  and  N2  as  necessary 

b  -  change  'DIMENSIONS'  according  to  new  values  of   N^  and  N2 
c  -  define  distribution  of  Mach  lines  in  the  Prandtl-ileyer  fan  as  required 
(program  lines  126-131) 

4.  Execution   commands: 

After   copying   program   into   USER'S   FILE: 

WATFIV      AXSYM   *    (XTYPE 
The   program   will    run   on   user's    terminal   under   WATFIV.      A  soft   copy   of 
program   listing   and   output   listing   will    be   stored    in   user's    disk   named 
AXSYM   LISTING. 
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5.  Hard   Copy 

PRINT   AXSYM  LISTING. 

6.  Program  Outputs. 

All  necessary  outputs  are  automatically  listed  by  the  program. 

Figure  (20)  and  Figure  (21)  show  the  resulting  mesh  of  characteristics 
calculated  for  an  altitude  of  200  km  for  110=4  and  M0=2  respectively. 

In  these  figures  we  also  show  some  isotherms  and  the  limit  where  the 
breakdown  parameter  equals  0.05.   These  lines  are  plotted  (manually)  using 
interpolation  procedures.   Data  along  the  breakdown  line  is  input  data  for 
the  molecular  flow. 
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APPENDIX  B   SIMUL  PROGRAM 

B.l   DATA  ORGANIZATION 

Because  of  the  large  number  of  molecules,  cells,  regions  and  sectors  in 
the  simulation  and  the  large  number  of  data  related  to  each  molecule,  each 
cell  and  region  to  be  stored,  special  precautions  should  be  taken  in  order  not 
to  overflow  the  available  computer  memory. 

The  following  data  organization  was  used  in  SII1UL.   Figure  22  shows  tne 
geometry  of  one  sector. 
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Figure  22.   Definition  of  sector  geometry  and  cell  volume. 

cell  volume  =  v(  I)  =  R(  I)  *  DDALFA  x  DR  x  RSI  *  DFI, 

v(I)    R(I)  •  RSI(I) 
v(l)    R(l)  .RSI 

v(l)  is  the  volume  of  smallest  cell  in  a  region. 
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1 .  Tables  of  Molecules  and  Their  Parameters 
P1(L,N1M)  -  Light  molecules  of  the  jet 
P2(L,N2M)   -   Heavy  molecules  of  the  jet 

P3(L,N3M)   -  Ambient  gas  molecules  (not  used  in  the  present  program) 
N1M,  N2M,  N3M  -  maximum  number  of  molecules  in  simulation.   Number  of 

active  molecules  may  be  smaller  or  equal  to  (N<I>tl). 
L=l,2,3   -   cartesian  components  of  velocity  vx,  vy,  vz  [m/s] 
L=4   -   radial  coordinate  [meters] 

if   r=-99  the  particular  molecule  is  inactive 
L=5  -  angular  coordinate  [radians]. 
This  table  is  generated  each  time  the  simulation  is  initiated  for  a  region. 
That  means,  the  same  group  of  molecules  (as  stored  in  the  computer  memory)  is 
used  to  simulate  the  flow  in  all  regions  in  the  computation  domain. 
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2.      C(M,I,j)    Table    of    Cells    (in   a  Region)-Real   Data 

I  =    1,10       -     radial  index  of   the   cell  in  a  region 

j    =    1,10        -     angular    index   of   the   cell   in  a  region 

M  =    19  -     radial   coordinate    of   cell   center 

M  =   20  -     angular    coordinate    of   cell   center 

M  =    1,9  -     time   parameter    for    collisions    of    different    species    in   a 

cell 
M  =    10-18      -     maximum  relative    velocity  expected   for    collisions    of 
different    species    in  a  cell 


64 


3.   IP(N1A+N2A+N3A)  Table  of  the  addresses  of  the  active  molecules 

arranged  in  order  of  their  species  and  in  the  order  of  their  cells 
IC(N,I,j)   -   table  of  cells  (in  a  region)  integer  data 
I  =  1,10  -  radial  index 
j  =  1,10   -   angular  index 

N  =  1  -  number  of  molecules  (spec  1) 
N  =  2  -  number  of  molecules  (spec  2) 
N  =  3  -  number  of  molecules  (spec  3) 
M  =  4  -   (starting  address  -  1)  of  molecules  as  ordered  in  (IP) 
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4.   Reg(N,kR,kS)  Data  table  for  a  specific  region  (real) 

kR  =  1,10  -  index  of  a  region  in  a  sector 

kS  =  1,20   -  index  of  the  sector 

N  =  1  -  Dfl  =  differential  angle  (axisymmetric) 

DF1  is  a  weighting  factor 

N  =  2  -  DN1  =  number  density  (species  1) 

N  =  3  -  DN2  =  number  density  (species  1) 

N  =  4  -  VOL  1  =  actual  volume  of  smallest  cell  in  kR 

N  =  5  -  AREA1  INPUT  area  of  smallest  cell  in  kR 
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5.   Region  Geometry  and  Input  Flux 

R(10)   -   polar  radius  of  a  cell  in  a  region 
input  area  of  a  cell 


A(10)   - 
VOL(IO)   - 


input  area  of  smallest  cell 
volume  of  a  cell 


volume  of  smallest  cell 
Ml (10)   -  initial  number  of  molecules  in  cells 

(equal  number  of  molecules  of  either  jet  species) 
Fl(10)  ,F2(10)   -   input  flux  through  high  pressure  starting  line 

(spec  1),  (spec  2) 
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6.  Input  Flux  From  High  Pressure  Starting  Line 
FWPl(N,I,j) 

FWP2(N,I,j) 
|  species 

positive  (input  molecules) 
west 
flux 
I  =  1,10   -   number  of  the  cell  along  the  starting  line  in  a  region 

(j) 

j  =  1,10   -   number  of  the  region  along  the  starting  line 
N  =  1   -  molecular  flow  for  a  given  DF1 
N  =  2  -  mean  molecular  velocity  Vx 
N  =  3   -  mean  molecular  velocity  Vy 
N  =  4  -  mean  input  gas  temperature 

7 .  Output  Flux 

FNNl(4,NAD,kR)  FSN1 (4 ,NAD,kR) 


FNN2(4,NAD,kR) 
ft 
| negative  (output) 

north 


FSN2(4,NAD,kR) 
I  negative  (output) 
south 


kR  -   number  of  region  in  the  sector 
NAD  -   angular  location 
The  first  index  include  the  same  parameters  as  (FWP1) 
FEN1(4,NRD)  FWN1(4,NRD) 


FEN2(4,NRD) 
ti 
| negative 


eas  t 


FWN2(4,NRD) 
ft 
negative 

wes  t 
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NRD  +  radial  location 

FEN1 ,FEN2 ,FWN1 ,FWN2  are  necessary  for  iterations  within  one  sector. 

8.   Sampling  of  Output  Flux  from  a  Sector 

After  averaging  they  are  transferred  to  F0E1 ,  F0E2 ,  F0W1 ,  F0W2. 

F0Wl(4,kC,kR,kS) 

output   flow   to    the  west 
F0W2(4,kC,kR,kS) 

F0El(4,kC,kR,kS) 

output  flow  to  the  east 
F0E2(4,kC,kR,kS) 

kC  -  radial  location  (cell)  of  flow  in  a  region 

kR  -   index  of  a  region  in  a  sector 

kS  -   index  of  the  sector 

The  first  index  (4  parameters)  include  the  same  parameters  as  FEP1(   ) 
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B.2   MOLECULAR  SIMULATION  FOR  A  GIVEN  REGION 

After  we  define  the  geometry  of  the  whole  sector  resulting  from  the 
region  geometry  and  cell  geometry,  we  may  start  with  the  molecular  simulation. 
This  includes : 

-  initial  setting  of  molecules  in  cells 

-  molecules  are  moved  according  to  the  time  increment  DTK 

-  new  molecules  are  generated  according  to  input  (or  output)  flows 

-  collisions  calculations 

-  integration  of  flow  parameters  for  average  parameters  calculations 

-  repetition  of  the  whole  procedure  as  long  as  necessary  to  obtain 
reasonable  statistical  averaging 

-  calculate  averages  and  flow  weighting. 

These  routines  are  the  core  of  the  program  and  must  be  repeated  for  all 
regions  in  a  sector  and  (or  all  sectors  in  the  domain  where  the 
collisions  are  significant). 

In  the  following  sections  we  bring  a  detailed  description  of  this  part  of 
the  program. 

1 .   Initial  Setting  of  Molecules  in  Cells 

a.  The  initial  number  of  available  molecules  in  a  simulation  is 
larger  than  the  number  of  active  molecules  (Pl(data,  number  of 
molecules),  P2(data,  number  of  molecules)  are  the  vectors  used  for 
species  1  and  2) 

b.  an  inactive  molecule  is  defined  as 

[P1(4,N)  or  P2(4,N)]  =  -99 
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c.  calculation  of  number  of  molecules  to  be  set  in  each  cell 

d.  calculation  of  cell  coordinates 

e.  deactivation  of  all  available  molecules  in  simulation 

f.  definition  of  molecules  coordinates 
P1(4,N),  P2(4,M)   are  polar  radiuses 

P1(5,N),  P2(5,M)   are  angular  coordinate  in  radians 

All  molecules  in  a  cell  are  set  at  random  locations  within  the 
cell. 

g.  definition  of  molecular  velocities 
P1(1,N),  P2(1,M)  velocity  in  X  direction 
P1(2,N),  P2(2,M)  velocity  in  Y  direction 
P1(3,N),  P2(3,H)  velocity  in  Z  direction 

Thermal  velocities  are  random  function  of  temperatures  and  are 
added  to  the  mean  velocity  as  defined  at  initial  boundary  ALFA0  of 
the  region. 
As  the  thermal  velocity  has  a  Boltzmann  distribution  the  thermal  velocity 
setting  is  based  on  rejection-acceptance  methods  (for  more  details  see  Bird 
[4]   Appendix  D) . 

h.   reset  collision  timers  and  relative  velocity 


71 


i.   reset  general  time  counter:   Time  =  0 
3 .   The  Simulation 

a.  Move  all  molecules  according  to  their  velocity  (Vx,  Vy)  and  find 
their  new  coordinates 

Note  -  A  routine  designed  to  calculate  the  collisions  with  the 
wall  was  included  in  the  program;  if  the  region  (or  sector)  is 
bounded  by  the  solid  wall,  the  collision  is  calculated  -  resulting 
new  velocities  and  directions  and  counted  for  wall  flux 
calculations.   If  the  program  is  stopped  at  an  angle  where  the 
flow  becomes  collisionless ,  this  routine  becomes  irrelevant  and 
other  type  of  calculations  should  be  designed. 

b.  Output  flow  counting:   all  molecules  that  leave  the  region  are 
counted  and  stored  in  specific  vectors  which  are  used  as 

inputs  to  other  regions.   The  output  vectors  are  (X  represent  1  or 

2  for  the  two  species  in  the  program  (FSNX)) 

FSNX(l,j,kR)   >   "south"  boundary 

FNNX(l,j,kR)  ■>      "north"  boundary 

FENX(1,I)      >  "east"   boundary 

FWNX(1,I)      •»■   "west"  boundary 

I  represents  the  radial  location  index  of  the  cell 

j  represents  the  angular  location  index 

kR  is  the  index  of  the  region  within  a  sector 
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Note  -  all  molecules  that  move  to 

(j=NAD+l)(I=NRD+l)  are  placed  in  FENX(1,NRD) 

(j=  -1)(X=NRD+1)  are  placed  in  FWNX(1,NRD) 

(j=NAD+l)(I=  -1)  are  placed  in  FSNX(l,l,kR) 

(j=  -1)(I=  -1)  are  placed  in  FWNX(l,l,kR) 
This  was  done  only  for  simplification  reasons. 

Generation  of  new  molecules  due  to  input  flows.   Through  the  four 
sides  E,N,W,S,  of  a  region,  molecules  are  allowed  to  enter  the 
region  according  with  the  flows  coming  from  the  neighbouring 
regions : 
for  "W"  side  of  the  region  in  sector  1 


for  other  sectors 

for  "E"  side  of  any  region 

for  "N"  side  of  any  region 

for  "S"  side  of  the  region 


FWPX(P,I,kR) 
FOEX(P,I,kR,kS-l) 
FOWX(P,I,kR,kS+l) 
FSNX(P,j ,kR+l) 


FNNX(P,j , kR- 1 ) 
The  first  parameter  of  all  these  arrays  represent: 
P  =  1  -  number  flux  (real  number) 
P  =  2  -  velocity  component  -  (Vx) 
P  =  3  -  velocity  component  -  (Vy) 
P  =  4  -  gas  temperature 
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Note  1  -  because  every  region  has  a  different  size  of  angle  DFI,  all 

fluxes  have  to  be  adjusted  accordingly. 

Note  2  -  input  fluxes  are  calculated,  adjusted  and  stored  as  real 

numbers.   The  number  of  input  molecules  are  by  definition  integers. 

In  order  not  to  "loose"  molecules,  the  number  of  input  molecules  is 

increased  by  1  on  a  random  basis.   (The  average  of  many  runs  will 

result  in  the  accurate  average  input  flow.) 
New  molecules  are  set  at  random  locations  on  the  boundary  of  the  specific 
cells  and  at  random  time  within  DTK.   Then  each  molecule  is  allowed  to  enter 
the  region  according  to  its  initial  coordinate  and  velocity.   At  the  end  of 
the  time  interval  the  new  location  and  velocity  is  stored  in  molecule  array 
PI  or  P2 .   If  DTM  is  chosen  to  be  too  large  and  cell  size  is  small  (total 
region  size  too  small)  some  molecules  may  cross  the  region  and  will  not  be 
counted  in  the  simulation  of  the  specific  region.   In  order  not  to  "loose" 
molecules : 

(a)  DTM  should  be  decreased 

(b)  count  these  molecules  as  output  fluxes  from  the  specific 
region.   (This  is  recommended  only  if  there  is  no  other 
choice. ) 

Note  -  Because  the  arrays  of  input  flows  store  only  averaged  data  for 

the  molecules,  the  thermal  velocity  of  each  new  molecule  is  calculated 

according  to  the  Boltzraann  distribution  as  a  function  of  the  averaged 

temperature. 

d.   Rearrangement  of  molecules  in  cells.   Before  collisions  are 

calculated  all  simulated  molecules  which  have  been  let  to  move  and 
generated  have  to  be  rearranged  and  recounted  for  each  cell. 
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The  array  IC(k,I,j)  contains  integer  data  for  each  cell 

(I.j) 

k  «■  1  -  is  the  number  of  molecules  of  species  1 

k  =  2  -  is  the  number  of  molecules  of  species  2 

k  ■  3  -  is  the  number  of  molecules  of  species  3  (not  used) 

k  ■  4  -  is  the  (address-1)  of  the  first  molecule  in  the  cell 

related  to  the  vector  IP(M) 

Vector  IP(M)  contains  the  list  of  the  simulated  molecules  arranged 

in  the  order  of  species  in  cells  and  cells  respectively.   The 

following  is  a  graphic  description  of  IP(M) 


<£.<£? Aj   7".cW  7~£ 


[RDDR£SS  -  -f ) 
of  Me  pA.it 

stored  <n  -rcfe,/,?) 


Figure   23.      Vector    IP(H) 
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B.3   SIMUL  Program  Flowchart 

A  simplified  flowchart  for  the  Monte  Carlo  simulation  of  the  molecular 
flow  is  shown  in  Figure  24.   The  program  is  designed  to  solve  the  ring 
axisymmetric  jet  flow,  however,  minor  changes  may  be  done  to  enable  a 
different  geometry. 
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Figure   24.      SIMUL   program  flowchart. 
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B.4   SIMUL  Program  Listing 

C      PROGRAM  SIMUL  SIM00010 

C  THIS  PROGRAM  IS  DESIGNED  TO  CALCULATE  THE  MOLECULAR  FLOW  OUTSIDE  THE  SIM00020 
C  CONTINUUM  REGION  FLOW  OF  A  HIGHLY  UNDEREXPANDED  AXISYMMETRIC  RING  JET.SIM00030 

C  RESULTS  FOR  THE  CONTINUUM  FLOW  MAY  BE  OBTAINED  FROM  'AXSYM'  PROGRAM  SIM00040 

C  WHICH  GIVES  THE  'METHOD  OF  CHARACTERISTICS'  ISENTROPIC  SOLUTION.  SIM00050 

C  THE  BOUNDARY  BETWEEN  THE  CONTINUUM  AND  MOLECULAR  FLOW  IS  DEFINED  BY  SIM00060 

C  THE  BREAKDOWN  PARAMETER  'P'  AS  PROPOSED  BY  G.  A.  BIRD.  SIM00070 

C  THE  'MOLECULAR  DOMAIN'  IS  DIVIDED  INTO  POLAR  DIVISIONS  MAKING  A  SET  SIM00080 

C  OF  SECTORS.  EACH  SECTOR  IS  SUDIVIDED  INTO  10  REGIONS  SIM00090 

C  NO  APRIORY  INFORMATION  ABOUT  THE  GEOMETRY  OF  THE  DIFFERENT  SIM00100 

C  REGIONS  IS  AVAILABLE,  THEREFORE,  MANUAL  INTERVENTION  MAY  BE  SIM00110 

C  REQUIRED  WHEN  MOVING  FROM  ONE  SECTOR  TO  THE  OTHER.  AN  ACCEPTABLE  SIM00120 

C  GEOMETRY  WILL  RESULT  A  REASONABLE  NUMBER  OF  SIMULATED  MOLECULES.  SIM00130 

C  EACH  REGION  IS  SUBDIVIDED  INTO  A  NUMBER  OF  CELLS  WITH  A  GEOMETRY  SIM0U140 

C  DEFINED  BY  A  POLAR  MESH.  A  NUMBER  OF  MOLECULES  IS  SET  IN  EACH  CELL  SIM00150 

C  PROPORTIONAL  TO  THE  CELL  VOLUME.  THE  BOUNDARY  CONDITIONS  FOR  EACH  SIM00160 

C  REGION  REQUIRE  INPUT  AND  OUTPUT  FLUX  OF  MOLECULES.  NO  APRIORY  SIM00170 

C  INFORMATION  ON  THE  FLUX  IS  AVAILABLE.  IT  WILL  BE  CALCULATED  IN  SIM00180 

C  AN  ITERATIVE  MODE.  SIM00190 

C  IF  THE  BOUNDARY  CONDITIONS  FOR  ALL  CELLS  ARE  CONSTANT  THE  NUMBER  SIM00200 

C  FLUX  INTO  EACH  CELL  IS  PROPORTIONAL  TO  CELL  WALL  AREA.  SIM00210 

C  TO  DECREASE  THE  ERROR  WHEN  INTRODUCING  MOLECULES  -(INTEGER  NUMBER)-  SIM00220 

C  DUE  TO  THE  INPUT  FLUX  -(REAL  NUMBER),  AN  ADDITIONAL  MOLECULE  IS  SIM00230 

C  GENERATED  ON  THE  BASIS  OF  RANDOM  NUMBERS  SUCH  THAT  THE  AVERAGE  OF  A  SIM00240 

C  LARGE  NUMBER  OF  SAMPLINGS  WILL  EQUAL  THE  REAL  INPUT  FLUX.  SIM00250 

Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxkxxxxxxxxxxxxxxxxxxxx  SI MO 0  26  0 

C         THE  MONTE  CARLO  SIMULATION  SIM00270 

C  THE  JET  GAS  IS  COMPOSED  OF  TWO  SPECIES  OF  MOLECULES  SIM00280 

C  AMBIENT  GAS  IS  REGARDED  AS  ONE  SPECIES  SIM00290 

C  THE  MOLECULAR  MODEL  IS  -  'HARD  SPHERE  MOLECULE'  SIM00300 

C  NETWORK  DEFINITION  SIM00310 

C     THE  MAXIMUM  RADIUS  (POLAR)  OF  THE  DOMAIN  IS  ASSUMED  TO  BE  RP=15  M.  SIM00320 

C     THE  ANGLE  OF  THE  BOUNDARY  BETWEEN  CONTINUUM  AND  MOLECULAR  FLOW  SIM00330 

C     AND  THE  SOLID  WALL  IS  SIM00340 

C                ALFA   (CALCULATED  IN  PROGRAM  'AXSYM')  SIM00350 

C     DEFINE  A  POLAR  SECTOR  WITH  A  RADIUS  ' RP '  AND  AN  ANGLE  OF  SIM00360 

C                DALFA  =5XMFP*NAD/RP  SIM00370 

C     THIS  SECTOR  IS  SUBDIVIDED  INTO  A  NUMBER  OF  RADIAL  DIVISIONS  SIM00380 

C     MAKING  »N'  REGIONS  FOR  SIMULATION  CALCULATIONS  FOR  EACH  'DALFA'.  SIM00390 

C     EACH  REGION  IS  DEVIDED  INTO  SIM00400 

C         NAD  -ANGULAR  DIVISIONS  (15)  WITH  AN  ANGLE  OF  DDALFA=5*MFP/RP  SIM00410 

C         NRD  -RADIAL   DIVISIONS  (10)  SIM00420 

C             MAKING  NADXNRD  CELLS  SIM00430 

C     THE  SMALLEST  CELL  CONTAINS  »MIN»  MOLECULES  SIM00440 

C                    MIN  =  15  SIM00450 

C     TOTAL  NUMBER  OF  ACTIVE  MOLECULES  IN  A  REGION  IS  LIMITED  TO  6000  SIM00460 

C     (3000  MOLECULES  OF  EACH  SPECIES.)  SIM00470 

C  THE  WIDTH  OF  A  CELL  DEFINED  BY  ANGLE  'DFI»  MAY  NOW  BE  EVALUATED.  SIM00480 

C  (MIN/NUMBER  DENSITY  =  ACTUAL  CELL  VOLUME)  SIM00490 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM0  0500 

C        DEFINE  PARAMETERS BLOCK  1SIM00510 

COMMON  IX  SIM00520 

DIMENSION  P1(5,3000),P2(5,3000),IP(6000),C(20,10,15),  SIM00530 

*IC(4,10,15),REG(5,10,20)  SIM00540 

C      DIMENSION  DFI(10),DN(10)  SIM00550 

DIMENSION  R(10),A(10),VOL(10),M1(10),F1(10),F2(10)  SIM00560 

DIMENSION  SS(2,10,10)  SIM00570 

DIMENSION  FEN1(<4,10),FWN1(<4,10),FNN1(4,15,10),FSN1(4,15,10)  SIM0  058  0 

DIMENSION  FEN2(4,10),FWN2(4,10),FNN2(4,15,10),FSN2(4,15,10)  SIM0  0590 

DIMENSION  FWPl(4,10,10),FWP2C<t,10,10)  SIM00600 

DIMENSION  FOE1(4,10,10,20),FOW1(4,10,10,20)  SIM00610 

DIMENSION  FOE2(4,10,10,20),FOW2(4,10,10,20)  SIM00620 

DIMENSION  NM(3),VRC(3),TOC(3,3),SPEC(3,5),NCOL(10,15)  SIM00630 

C  P1,P2,P3  CONTAIN  INFORMATION  ON  UP  TO  M  SIMULATED  MOLECULES  SIM00640 

C    P1(1,N),P1(2,N),P1(3,N)  ARE  U,V,W  VELOCITY  COMPONENTS  (CARTESIAN)  SIM00650 

C    P1(4,N),P1(5,N)  ARE  R  AND  TETA  COORDINATES  (POLAR)  SIM00660 

C  SIM00670 

C  IP(M)  ARE  THE  M  MOLECULES  ARRANGED  IN  THE  ORDER  OF  THEIR  CELLS  SIM00680 

C  C(20,I,J)  CONTAINS  INFORMATION  ON  UP  TO  IXJ  CELLS  SIM00690 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM00700 

C  SIM00710 

Z    ALFA  IS  THE  ANGLE  WHERE  THE  BREAKDOWN  PARAMETER  EQUALS  .05    .  SIM00720 
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C  FPM  IS  THE  MEAN  FREE  PATH  AT  ALFA.  SIM00730 
JCXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIMO 074  0 

C  SET  GENERAL  CONSTANTS BLOCK  2SIM00750 

i      IX     =  529814367  SIM00760 

PI      =  3.141593  SIM00770 

BOLTZ   =  1.38044E-23  SIM00780 

PO      =  101325.  SIM00790 

TO      =  273.  SIM00800 

AVOG    =  2.68699E+25  SIM00810 

RG      =  8314.  SIM00820 
CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxsiM0  0830 

C  SET  PROGRAM  CONSTANTS BLOCK  3SIM00840 

ISPEC=2  SIM00850 

Rl  =  2.5  SIM00860 

DR=.15  SIM00870 

RP=15.  SIM00880 

C  NUMBER  OF  SIMULATED  MOLECULES  IN  SMALLEST  CELL  (DIFFERENT  SPECIES)     SIM00890 

NMOL1=3000  SIM00900 

NMOL2=3000  SIM00910 

NMOL3=0  SIM00920 

MIN  =  15  SIM00930 

N1M  =  15  SIM00940 

N2M  =  15  SIM00950 

C  NUMBER  OF  DIVISIONS  IN  A  SIMULATED  REGION  SIM00960 

NAD  =  15  SIM00970 

NRDS=10  SIM00980 

NRD=10  SIM00990 

NRD1=  3  SIM01000 

NIS  =2  SIM01010 

C  Rl  IS  THE  RADIUS  OF  THE  CYLINDER  (WALL)  SIM01020 

C  R(I)  IS  THE  RADIAL  COORDINATE  OF  CELL(I)  SIM01030 

C  VOL(I)  IS  THE  RATIO  BETWEEN  VOLUMES  OF  CELL(I)  AND  CELL(l)  SIM01040 

C  A(I)  IS  THE  RATIO  BETWEEN  INPUT  FLOW  AREAS  OF  CELL(I)  AND  CELL(l)      SIM01050 

C  DR  IS  THE  RADIAL  SIZE  OF  A  CELL  (RADIAL  INCREMENT)  SIM01060 

C  INPUT  THE  FOLLOWING  AS  'DATA1  OR  'READ*  STATEMENTS  xxxxxxxxxxxxxxxx    SIM01070 

TETA  =  1.2  SIM01080 

C  TETA  IS  THE  FLOW  DIRECTION  (RADIANS)  ON  THE  BREAKDOWN  LINE  AS  FOUND    SIM01090 

C  FROM  THE  AXSYM  PROGRAM.  VO  IS  THE  FLOW  VELOCITY  (M/SEC)  SIM01100 

VO    =  2100  SIM01110 

VOX   =  VOXCOS(TETA)  SIM01120 

VOY   =  VOXSIN(TETA)  SIM01130 

TWALL=300.  SIM01140 

DN01=1.1E21  SIM01150 

DN02=1.1E21  SIM01160 

DTM=1.5E-6  SIM01170 

C  SET  DATA  FOR  THE  DIFFERENT  SPECIES  SIM01180 

C    MOLECULAR  MASS  SIM01190 

SPEC(1,1)=4./AV0G  SIM01200 

SPEC(2,1)=40./AV0G  SIM01210 

SPEC(3,1)=29./AV0G  SIM01220 

C    MOLECULAR  DIAMETER  SIM01230 

SPEC(1,2)=2.19E-10  SIM01240 

SPEC(2,2)=4.00E-10  SIM01250 

C      SPEC(3,2)=  SIM01260 

C     SPEC(1,3)=  SIM01270 

C      SPEC(2,3)=      ADDITIONAL  DATA  IF  REQUIRED  SIM01280 

C      SPEC(3,3)=  SIM01290 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM0130  0 

C  THERMAL  VELOCITIES  AT  WALL  TEMPERATURE  SIM01310 

VWM1=SQRT(2.*B0LTZXTWALL/SPEC(1,1))  SIM0132  0 

VWM2  =  SQRT(2.XB0LTZXTWALL/SPEC(2,D)  SIM01330 

C  SIM01340 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM01350 

C  INITIALISATION. .RESET  ALL  SAMPLING  VARIABLES  SIM01360 

DO  80  1  =  1, NRD  SIM01370 

DO  80  JR=1,10  SIM01380 

DO  80  JS=1,20  SIM01390 

DO  80  KPAR=1,4  SIM01400 

F0E1(KPAR,I,JR,JS)=0.  SIM01410 

F0E2(KPAR,I, JR, JS)=0.  SIM01420 

F0W1(KPAR,I,JR, JS)=0.  SIM01430 

80  F0W2(KPAR,I,JR,JS)=0.  SIM01440 
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C  SET  INPUT  PARAMETERS BLOCK  3  SIM0146C 

ITER=1  SIM0147C 

C  RETURN  TO  3000  FOR  ADDITIONAL  ITERATIONS  SIM0148C 

3000  KS=1  SIM0149C 

ALFA=1.3  SIM0150C 

FPM=.001  SIM0151C 

TEMP=40.  SIM0152( 

VTER1=SQRT(2.*B0LTZ*TEMP/SPEC(1,1))  SIM0153C 

VTER2=SQRT(2.*B0LTZXTEMP/SPEC(2,1))  SIM0154C 

IF  (ITER.GT.l)GO  TO  2005  SIM0155C 

DO  2001  1=1,10  SIM0156( 

REG(2,I,1)=DN01  SIM0157( 

REG(3,I,1)=DN02  SIM0158( 

DO  2001  J=l,10  SIM0159( 

FWP1(2,I,J)=V0X  SIM0160( 

FWP2(2,I,J)=V0X  SIM016K 

FWP1C3,I,J)=V0Y  SIM01621 

FWP2C3,I,J)=V0Y  SIM0163( 

FWP1(4,I,J)=TEMP  SIM0164( 

2001  FWP2(4,I,J)=TEMP  SIM0165( 

2005  CONTINUE  SIM01661 

C  SIM0167I 

C3€XXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX9(XXXXXXXXXSIM0168( 

C  SIM0169I 

CXXXXXXXXXXXXXX«XXXXXXXXXXKXXXXXXXXXXXX»XXXXXXX»XXXXXXXXXXXXKXXXXXXKXXXXSIM017  0( 

C  DEFINE  SECTOR  GEOMETRY .  .  .BLOCK  <♦  SIM017K 

C  KR  IS  THE  INDEX  FOR  A  REGION  IN  ONE  SECTOR  SIM01721 

C  RETURN  TO  1000  FOR  NEXT  SECTOR  SIM0173I 

DALFA=0.  SIM01741 

1000  KR=1  SIM0175I 

C  RESET  SAMPLING  OF  FLOW  VARIABLES     (N  AND  S)  SIM0176I 

DO  81  1=1, NAD  SIM0177I 

DO  81  JR=1,10  SIM0178I 

DO  81  KPAR=1,4  SIM0179I 

FNN1(KPAR,I, JR)=0.  SIM0180I 

FNN2(KPAR,I, JR)=0.  SIM018K 

FSN1(KPAR,I,JR)=0.  SIM0182I 

81  FSN2(KPAR,I,JR)=0.  SIM0183I 
C  ALFA  IS  THE  STARTING  ANGLE  OF  THE  PRESENT  SECTOR.  PREVIOUS  DALFA  SIM0184I 
C  SUBTRACTED  FROM  PREVIOUS  ALFA  (PREVIOUS  SECTOR)  GIVES  THE  NEW  ALFA.    SIM0185I 

ALFA=ALFA-DALFA  SIM0186I 

DALFA  =  5.XFPM*FL0AT(NAD)/RP  SIM0187I 

C  DALFA  IS  THE  ANGLE  OF  THE  SECTOR  SIM0188( 

IFCDALFA.GT. ALFA/2. ) DALFA=ALFA/2 .  SIM0189I 

IF(DALFA.GT.(FPM*75/RP))DALFA=ALFA  SIM0190I 

IFCALFA.LE.O. )GO  TO  2000  SIM019U 

DDALFA  =  DALFA/FLOATCNAD)  SIM0192I 

C  DDALFA  IS  THE  ANGLE  OF  A  SINGLE  CELL  IN  THE  SECTOR                    SIM0193I 

ALFAJ=ALFA-DALFA/2.  SIM0194( 

C  SIM01951 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXSIM0196I 

C  DEFINE  REGION  KR  IN  SECTOR BLOCK  5  SIM0197I 

C  RETURN  TO  2000  FOR  THE  NEXT  REGION  KR  SIM0198( 

2000  CONTINUE  SIM0199( 

C  RESET  SAMPLING  OF  FLOW  VARIABLES  PER  REGION  SIM0200( 

DO  82  1=1, NRD  SIM020K 

DO  82  KPAR=1,4  SIM0202( 

FEN1(KPAR,I)=0.  SIM0203C 

FEN2(KPAR,I)=0.  SIM0204C 

FWN1(KPAR,I)=0.  SIM0205C 

82  FWN2(KPAR,I)=0.  SIM0206C 
MT=0  SIM0207C 
NRD=NRDS  SIM0208C 
IF(KR.LT.2)NRD=NRD1  SIM0209C 
DO  100  I  =  1,NRD  SIM0210C 

C  POLAR  RADIUS  MEASURED  FROM  THE  NOZZLE  LIP  SIM0211C 

R(I)  =(FL0AT(I)-.5)*DR  SIM0212C 

IF(KR.EQ.2)R(I)=R(I)+FL0AT(NRD1)*DR  SIM0213C 

IF(KR.GE.3)R(I)=R(I)+(FL0AT(NRD1+(KR-2)XNRD))*DR  SIM0214C 

C  RSI  IS  THE  RADIUS  MEASURED  FROM  THE  AXIS  OF  SYMMETRY                  SIM0215C 

RSI=R(I)+R1/SIN(ALFAJ)  SIM0216C 
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IF(I.EQ.1)RS1=RSI  SIM02170 

A(I)    =  RSI/RSI  SIM02180 

VOL(I)  =  R(I)XRSI/(R(1)XRS1)  SIM02190 

M1CI)    =  MINXVOL(I)  SIM02200 

MT  =  MT+M1CI)  SIM02210 

C  MK)  IS  THE  INITIAL  NUMBER  OF  MOLECULES  IN  EACH  CELL  SIM02220 

DN1=REG(2,KR,KS)  SIM02230 

DN2=REG(3,KR,KS)  SIM02240 

REG(l,KR,KS)=FLOAT(MIN)/(DNlxR(l)XRSlxDRxDDALFA)  SIM02250 

C  DFI  (REG(1,KR,KS))  IS  THE  SPHERICAL  ( AXISYMMETRIC)  ANGLE  SIM02260 

C  THIS  ANGLE  HAS  DIFFERENT  VALUES  FOR  EACH  KR,  THEREFORE  IT  IS  A         SIM02270 

C  WEIGHTING  FACTOR .+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+++++++++++++   SIM02280 

C   (FNN1,FNN2,FSN1,FSN2,FEN1,FEN2,FWN1,FWN2)  X    DFI  SIM02290 

C   CF0E1,F0E2,F0WI,F0W2)                      X    DFI  SIM02300 

DAI  =  RS1XDRXREG(1,KR,KS)  SIM02310 

C  DAI  IS  THE  INPUT  AREA  OF  CELL(l)  SIM02320 

FN1  =  V0XDA1XREG(2,KR,KS)XSIN(ALFA-TETA)  SIM02330 

FN2  =  V0XDA1XREG(3,KR,KS)XSIN(ALFA-TETA)  SIM02340 

C  TETA  IS  THE  ANGLE  BETWEEN  FLOW  DIRECTION  AND  THE  WALL  SIMQ2350 

FKI)  =  FNIXA(I)  SIM02360 

F2(I)  =  FN2XACI)  SIM02370 

C  SIM02380 

99  FWPK1,I,KR)=FKI)  SIM02390 

100  FWP2(1,I,KR)=F2(I)  SIM02400 

C SIM02410 

C  DEFINE  ACTUAL  VOLUME  AND  INPUT  AREA  OF  SMALLEST  CELL  IN  REGION  SIM02420 

REG(4,KR,KS)=R(1)XDDALFAXDRXRS1XREG(1,KR,KS)  SIM02430 

REG(5,KR,KS)=DA1  SIM02440 

C  FN1  IS  THE  NUMBER  FLUX  TO  THE  SMALLEST  CELL  (REAL  NUMBER)  PER  SECOND.  SIM02450 

C  INPUT  NUMBER  OF  MOLECULES  (INTEGER)  WILL  BE  INTEGRATED  TO  MAKE  AN      SIM02460 

C  AVERAGE  OF  FN1 .  SIM02470 

C  FKI)  IS  THE  INPUT  FLUX  TO  CELL(I)  SIM02480 

C   FOLLOWING  ARE  THE  REGION  BOUNDARIES  SIM02490 

102  RMIN=R(1)-.5XDR  SIM02500 
RMAX=RMIN+DRXFL0AT(NRD)  SIM02510 
TMAX=ALFA  •  SIM02520 
TMIN=ALFA-DALFA  SIM02530 
DO  9102  I  =  1,NRD  SIM02540 
WRITE  (6,9101)R(I),A(I),VOL(I),MKI),FKI)  SIM02550 

9101  FORMATC  » , 3F13 . 5, 110 , E15 . 5)  SIM02560 

9102  CONTINUE  SIM02570 
MTT=  MTXNAD  SIM02580 
WRITE(6,103)MT,MTT  SIM02590 

103  FORMAT  ('  INITIAL  NUMBER  OF  SIMULATED  MOLECULES  IS=  ',15,'  PER  SPESIM02600 
XCIES  PER  DDALFA,  TOTAL  NUMBER  IS',  15)  SIM02610 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM0  26  20 

C  DEFINE  INPUT  FLOWS  TO  KR  W,E,N,S BLOCK  6  SIM02630 

C  SIM02640 

C  SIM02650 

C  SIM02660 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM0267  0 

C  CALCULATE  INITIAL  NUMBER  OF  MOLECULES  IN  CELLS BLOCK  7  SIM02680 

C  SET  INITIAL  STATE  OF  GAS •  »   8  SIM02690 

C  SIM02700 

DO  150  I  =1,NRD  SIM02710 

DO  150  J=  1,NAD  SIM02720 

C  SIM02730 

C  SET  SIMULATED  MOLECULES  IN  THEIR  CELLS  SIM02740 

IC(1,I,J)  =  MKI)  SIM02750 

IC(2,I,J)  =  MKI)  SIM02760 

IC(3,I,J)  =  0  SIM02770 

C  SIM02780 

C  SET  CELL  COORDINATES  SIM02790 

C(19,I,J)  =  R(I)  SIM02800 

C(20,I,J)  =  ALFA-(FL0AT(J)-.5)XDDALFA  SIM02810 

150  CONTINUE  SIM02820 

Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx  SI MO 28 30 

C  SIM028^0 

C  DEACTIVATE  ALL  MOLECULES  SIM02850 

DO  170  N  =  KNMOL1  SIM02860 

PK4,N)  =  -99.  SIM02870 

170  P2(4,N)  =  -99.  SIM02880 
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C  SIM02890 

Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx  SIM02900 

C  SET  INITIAL  STATE  OF  THE  GAS  (LOCATION  AND  VELOCITY  OF  MOLECULES)  SIM02910 

C  SIM02920 

NADR1  =  0  SIM02930 

NADR2  =  0  SIM02940 

DO  200  I  =  1,NRD  SIM02950 

DO  200  J  =  1,NAD  SIM02960 

NM1  =  ICC1,I,J)  SIM02970 

C  SIM02980 

DO  205  N  =  1,NM1  SIM02990 

NADR1  =  NADR1  +  1  SIM03000 

CALL  RANDU(PP)  SIM03010 

P1(4,NADR1)  =  C(19,I,J)+DRX(PP-„5)  SIM03020 

CALL  RANDU(P)  SIM03030 

P1C5,NADR1)  =  C(20,I,J)+DDALFAX(P-.5)  SIM03040 

C  SIM03050 

DO  205  NV  =  1,3  SIM03060 

203  CALL  RANDU(P)  SIM03070 

V  =  -3.+6.XP  SIM03080 
B  =  EXP(-VXV)  SIM03090 
CALL  RANDU(P)  SIM03100 
IF(B.LT.P)  GO  TO  203  SIM03110 
P1(NV,NADR1)  =  P*SINCB)XVTER1  SIM03120 
IF(NV.EQ.1)P1CNV,NADR1)=P1(NV,NADR1)+V0X  SIM03130 
IF(NV.EQ.2)P1(NV,NADR1)=P1CNV,NADR1)+V0Y  SIM03140 

205  CONTINUE  SIM03150 

C  SIM03160 

C  REPEAT  PROCEDURE  FOR  SPECIES  2  SIM03170 

NM2  =IC(2,I,J)  SIM03180 

DO  210  N  =  1,NM2  SIM03190 

NADR2  =  NADR2+1  SIM03200 

CALL  RANDU(P)  SIM03210 

P2C4,NADR2)  =  C(19,I,J)+DRX(P-.5)  SIM03220 

CALL  RANDU(P)  SIM03230 

P2(5,NADR2)  =  C( 20 , I , J )+DDALFA*(P- . 5)  SIM03240 

C  SIM03250 

DO  210  NV  =  1,3  SIM03260 

207  CALL  RANDU(P)  SIM03270 

V  =  -3.+6.XP  SIM03230 
B  =  EXPC-VXV)  SIM03290 
CALL  RANDU(P)  SIM03300 
IF(B.LT.P)  GO  TO  207  SIM03310 
P2(NV,NADR2)  =  PXSINC B)XVTER2  SIM03320 
IF(NV.EQ.1)P2(NV,NADR2)=P2(NV,NADR2)+V0X  SIM03330 
IF(NV.EQ.2)P2(NV,NADR2)=P2(NV,NADR2)+V0Y  SIM03340 

210  CONTINUE  SIM03350 

C  SIM03360 

C  WHEN  NECESSARY,  REPEAT  PROCEDURE  FOR  SPECIES  3.  SIM03370 

200  CONTINUE  SIM03380 

C  SIM03390 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM0  3400 

C  DEFINE  HERE  ALL  COLLISION  PARAMETERS  TO  BE  INCLUDED  IN  SIMULATION  SIM03410 

C  BLOCK  9  SIM03420 

C  RESET  COLLISION  TIMERS  SIM03430 

DO  19  1  =  1,  NRD  SIM03440 

DO  19  J=1,NAD  SIM03450 

DO  19  L=l,3  SIM03460 

DO  19M=1,3  SIM03470 

KT=3X(L-l)+M+9  SIM03480 

C(KT,I,J)=0.  SIM03490 

C  SET  EXPECTED  MAXIMUM  RELATIVE  VELOCITY  IN  COLLISIONS  SIM03500 

KV=3X(L-1)+M  SIM03510 

EM1=SPECCL,1)  SIM03520 

EM2=SPEC(M,1)  SIM03530 

EMR=EM1XEM2/(EM1+EM2)  SIM03540 

C  SIM03550 

IFCKS.EQ.DGO    TO    17  SIM03560 

TEM=F0E1(4,I,KR,KS)  SIM03570 

GO  TO  18  SIM03580 

17  TEM=FWP1(4,I,KR)  SIM03590 

C  SIM03600 


18  RV=2./SQRTCPIX2.XB0LTZ*TEM/EMR)  SIM03610 

19  C(KV,I,J)=2.XRV  SIM03620 
C  SIM03630 
C  MAXIMUM  RELATIVE  VELOCITY  WILL  BE  RESET  IF  FASTER  ENCOUNTERS  OCCUR  SIM03640 
C  SIM03650 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM0  366  0 

TIME=0.  SIM03670 

C  LOOP  OVER  TIME  INTERVALS BLOCK  10  SIM03680 

DO  6000  JDTM  =  1,NIS  SIM03690 

C  SIM03700 

C  SIM03710 

TIME=TIME+DTM  SIM03720 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM037  30 

C  SIM03740 

C  MOVE  ALL  MOLECULES  BLOCK  12  SIM03750 

C  MOVE  MOLECULES  OF  SPECIES  1  SIM03760 

C  SIM03770 

DO  310  II  =  1,NM0L1  SIM03780 

C  SKIP  INACTIVE  MOLECULES 13  SIMJ53790 

IF  (Pl(«+,Il).EQ.-99.)  GO  TO  310  SIM03800 

VX  =  P1(1,I1)  SIM03810 

VY  =  P1(2,I1)  SIM03820 

RX  =  P1(<4,I1)  SIM03830 

T  =  P1(5,I1)  SIM03840 

TOLD=T  SIM03850 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM0386  0 

C  FIND  NEW  COORDINATES BLOCK  14  SIM03870 

XX=RXXC0S(T)+VXXDTM  SIM03880 

YY=RXXSIN(T)+VYXDTM  SIM03890 

RNEW=SQRT(XXXX2+YYX*2)  SIM03900 

T=ATAN(YY/XX)  SIM03910 
Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx  SIM03920 

C  FOR  LAST  SECTOR  FIND  COLLISIONS  WITH  THE  WALL  AND  SAMPLE  THEM  SIM03930 

C  BLOCK  15,16   SIM03940 

IF(ALFA.GT.DALFA)GO  TO  301  SIM03950 

IF(T.GT.0.)GO  TO  301  SIM03960 

DTR=DTMXT/(T-TOLD)  SIM03970 

C  DTR  IS  THE  TIME  REMAINING  AFTER  A  MOLECULE  STRIKES  THE  WALL    -  SIM03980 

IF(DTR.LT.1E-10)DTR=1E-10  SIM0399  0 

RW=RX+(RNEW-RX)XDTR/DTM  SIM04000 

IF(RW.LT.RMIN)RW=RMIN+DTRX.001  SIM04010 

IF(RW.GT.RMAX)RW=RMAX-DTRX.0  01  SIM04020 

LOC=RW/DR+l.  SIM04030 

C  SIM04040 

C  COUNT  COLLISIONS  WITH  THE  WALL(MUST  BE  WAIGHTED  =XDFI)  SIM04050 

SSC1,L0C, JSAMP)=SSC1,L0C,JSAMP)+1  SIM0406  0 

C  SIM04070 

C  SET  VELOCITY  AND  LOCATION  AFTER  A  MOLECULE  STRIKES  THE  WALL  SIM04080 

302  CALL  RANDU(P)  SIM04090 

B=VWM1XSQRT(-AL0G(P))  SIM04100 

CALL  RANDU(P)  SIM04110 

BB=2.XPIXP  SIM04120 

P1(1,NADR1)=BXC0S(BB)  SIM04130 

VX=BXC0S(BB)  SIM04140 

P1(2,NADR1)=BXSIN(BB)  SIM04150 

VY  =  BXSIN(BB)  SIM04160 

CALL  RANDU(P)  SIM04170 

P1(3,NADR1)=VWM1XSQRT(-AL0G(P))  SIM0418  0 

XX=RXXC0SCT)+VXXDTR  SIM04190 

YY=RXXSIN(T)+VYXDTR  SIM04200 

RNEW=SQRT(XXXX2+YYXX2)  SIM04210 

T  =  ATAN(YY/XX)  SIM04220 

C  SIM04230 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM0  4240 

C  DEACTIVATE  ALL  MOLECULES  THAT  MOVED  OUT BLOCK  17  SIM04250 

301  IF(RNEW.LT.RMAX.AND.RNEW.GT.RMIN.AND.T.LT.TMAX.AND.T.GT.TMIN)GO  T0SIM0426  0 

X303  SIM04270 

RI=(RNEW-RMIN)/DR+. 999999  SIM0428  0 

IR=RI  SIM04290 

TI=(TMAX-T)/DDALFA+. 999999  SIM0430  0 

IT=TI  SIM04310 

C  SIM04320 
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IFCIT.LE.0)GO  TO  304  SIM04330 

IF(IT.GT.NAD)GO  TO  305  SIM04340 

1       IF(IR.GT.NRD)GO  TO  306  SIM04350 

C  COUNT  S  DIRECTION,  SAMPLE 18   SIM04360 

FSN1C1,IT,KR)=FSN1C1,IT,KR)+1.  SIM0437  0 

GO  TO  309  SIM04380 

C  COUNT  W  DIRECTION,  SAMPLE. 18   SIM04390 

304  IF(IR.LE.0)IR=1  SIM04400 
IF(IR.GT.NRD)IR=NRD  SIM04410 
FWN1(1,IR)=FWN1C1,IR)+1.  SIM04420 
GO  TO  309  SIM04430 

C  COUNT  E  DIRECTION,  SAMPLE.  .  . 18   SIM04440 

305  IF(IR.LE.0)IR=1  SIM04450 
IF(IR.GT.NRD)IR=NRD  SIM04460 
FEN1(1,IR)=FEN1(1,IR)+1.  SIM0447  0 
GO  TO  309  SIM04480 

C  COUNT  N  DIRECTION,  SAMPLE .18   SIM04490 

306  FNN1(1,IT,KR)=FNN1(1,IT,KR)+1.  SIM04500 
C  SET  NEW  VALUES  IN  THE  MOLECULE  TABLE  SIM04510 
C  SIM04520 

309  RNEW=-99.  SIM04530 
303  P1(4,I1)=RNEW  SIM04540 

P1(5,I1)=T  SIM04550 

310  CONTINUE  SIM04560 
C  SIM04570 

C SIM04580 

C  REPEAT  PROCEDURE  FOR  SPECIES  2  MOLECULES BLOCKS  12  T018    SIM04590 

DO  320  I2=1,NM0L2  SIM04600 

C  SKIP  INACTIVE  MOLECULES  ,  SIM04610 

IF(P2C<t,I2).EQ.-99.)G0  TO  320  SIM04620 

VX  =  P2(1,I2)  SIM04630 

VY  =  P2(2,I2)  SIM04640 

RX  =  P2(4,I2)  SIM04650 

T  =  P2(5,I2)  SIM04660 

TOLD=T  SIMO4670 

C  SIM04680 

XX=RX*COSCT)+VXXDTM  SIM04690 

YY=RX*SIN(T)+VY*DTM  5IM0470C 

RNEW=SQRT(XXX*2+YY*X2)  SIM0471Q 

T=ATAN(YY/XX)  SIM04720 

C   COLLISIONS  WITH  THE  WALL  SIM0473Q 

Z. SIM0474G 

IF(ALFA.GT.DALFA)GO  TO  311  SIMO4750 

IFCT.GT.O. )GO  TO  311  SIMO4760 

DTR=DTM*T/(T-TOLD)  SIM04770 

C  DTR  IS  THE  TIME  REMAINING  AFTER  A  MOLECULE  STRIKES  THE  WALL  SIM04780 

IF(DTR.LT.1E-10)DTR=IE-10  SIM047  90 

RW=RX+(RNEW-RX)*DTR/DTM  SIM048  0  0 

IF(RW.LT.RMIN)RW=RW+DTRX.001  3IM04810 

IF(RW.GT.RMAX)RW=RW-DTRX.001  SIM04820 

LOC=RW/DR+l  SIM04830 

C  COUNT  COLLISIONS  WITH  THE  WALL  (MUST  BE  WEIGHTED  =*DFI)  SIM04840 

SS(2,L0C,JSAMP)=SS(2,L0C,JSAMP)+1  SIM04850 

C  SIM04860 

C  FIND  THE  NEW  COORDINATES  OF  THE  MOLECULE  SIM04S70 

312  CALL  RANDU(P)  SIMO4880 

B=VWM2*SQRT(-AL0G(P))  SIM04890 

CALL  RANDU(P)  SIM04900 

BB=2.*PI*P  SIM04910 

P2(1,NADR2)=B*C0S(BB)  SIM04920 

VX=B*COS(BB)  SIM04930 

P2(2,NADR2)=B*SIN(BB)  SIM04940 

VY=B*SIN(BB)  SIM04950 

CALL  RANDU(P)  SIM04960 

XX=RXXCOS(T)+VX*DTR  SIM04970 

YY=RX*SIN(T)+VY*DTR  SIM04980 

RNEW=SQRT(XX*X2+YY*X2)  SIM04990 

T=ATAN(YY/XX)  SIM05000 

C  SIM05010 

C  DEACTIVATE  MOLECULES  THAT  MOVED  OUT. COUNT  FOR  OUTPUT  FLUX  EVALUATION.  SIM05020 

311  IF(RNEW.LT.RMAX.AND.RNEW.GT.RMIN.AND.T.LT.TMAX.AND.T.GT.TMIN)GO  TOSIM0  50  30 
X313  SIM05040 
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RI  =  (RNEW-RMIN)/DR+. 999999  SIM05050 

IR  =  RI  SIM05060 

TI  =  <TMAX-T)/DDALFA+. 999999  SIM05070 

IT  =  TI  SIM05080 

IF  (IT.LE.O)GO  TO  314  SIM05090 

IF  (IT.GT.NAD)GO  TO  315  SIM05100 

IF(IR.GT.NRD)GO  TO  316  SIM05110 

C  COUNT  S  DIRECTION,  SAMPLE 18   SIM05120 

FSN2C1,IT,KR)=FSN2(1,IT,KR)+1.  SIM05130 

GO  TO  319  SIM05140 

C  COUNT  W  DIRECTION,  SAMPLE 18   SIM05150 

314  IF(IR.LE.O)IR=l  SIM05160 
IF(IR.GT.NRD)IR=NRD  SIM05170 
FWN2(1,IR)=FWN2(1,IR)+1.  SIM05180 
GO  TO  319  SIM05190 

C  COUNT  E  DIRECTION,  SAMPLE 18   SIM05200 

315  IF(IR.LE.0)IR=1  SIM05210 
IF(IR.GT.NRD)IR=NRD  SIM05220 
FEN2(1,IR)=FEN2C1,IR)+1.  SIM05230 
GO  TO  319  SIM05240 

C  COUNT  N  DIRECTION,  SAMPLE, 18   SIM05250 

316  FNN2(1,IT,KR)=FNN2C1,IT,KR)+1.  SIM05260 
C  SIM05270 

319  RNEW  =  -99.  SIM05280 
C 19   SIM0529  0 

313  P2(4,I2)=  RNEW  SIM05300 

P2(5,I2)=  T  SIM05310 

320  CONTINUE  SIM05320 
C  SIM05330 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM0  534  0 
C  NOW  NEW  MOLECULES  WILL  BE  INTRODUCED  SIM05350 
C  TOTAL  JET  FLUX  INTO  THE  REGION  WAS  DETERMINED  BY  FKI)  MOL./SEC.  SIM05360 
C  OR  FWP1(I,KR),FWP2(I,KR)  FOR  EACH  SPECIES  SIM05370 
C  NUMBER  OF  MOLECULES  TO  BE  ACTIVATED  PER  DTM  IS  FCI)XDTM  PER  CELL.  SIM05380 
C  SIM05390 
Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx   SIM05400 

DO  390  1=1, NRD  SIM05410 

C  INPUT  'W»  MOLECULES  SIM05420 

IFCKS.EQ.DGO  TO  322  SIM05430 

ANEW1=F0E1(1,I,KR,KS-1)  SIM05440 

SIM05450 
SIM05460 
SIM05470 
SIM05480 
SIM05490 
SIM05500 
SIM05510 
SIM05520 
SIM05530 
SIM05540 
SIM05550 
SIM05560 
SIM05570 
NEW  MOLS.  ARE  RANDOM  FUNCTIONS     SIM05580 

SIM05590 
SIM05600 

C  SIM05610 

DO  340  11=1, NEW1  SIM05620 

CALL  RANDU(P)  SIM05630 

ATIME=P*DTM  SIM05640 

CALL  RANDU(P)  SIM05650 

RSTART=(FLOAT(I-l)+P)*DR+RMIN  SIM0566  0 

VELX=FWP1(2,I,KR)  SIM05670 

VELY=FWP1(3,I,KR)  SIM05680 

TEM=FWP1(4,I,KR)  SIM05690 

C  THERMAL  VELOCITY  SIM05700 

VTER1=SQRT(2.XB0LTZXTEM/SPEC(1,1))  SIM05710 

C  TEMPERATURE  IS  TRANSFERRED  THROUGH  FWP1(4,I,KR)  SIM05720 

XX=RSTART*C0S(ALFA)+VELXXATIME  SIM057  30 

YY=RSTARTXSIN(ALFA)+VELYXATIME  SIM057  40 

RX=SQRT(XXXX2+YYXX2)  SIM05750 

T=ATAN(YY/XX)  SIM05760 
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ANEW2=F0E2(1,I,KR,KS-1) 

GO  TO  323 

322 

ANEW1=FWP1(1,I,KR)XDTM 
ANEW2=FWP2(1,I,KR)XDTM 

323 

NEW1=ANEW1 

NEW2=ANEW2 

REM1=ANEW1-NEW1 

CALL  RANDU(P) 

IF(REM1 .GT.P)NEW1=NEW1+1 

REM2=ANEW2-NEW2 

CALL  RANDU(P) 

IF(REM2.GT.P)NEW2=NEW2+1 

C 

ACTIVATE  NEW  INPUT  MOLECULES 

c 

TIME, LOCATION  AND  VELOCITY  COMP .  OF  N 

c 

SPE( 

:ies  l 

IF(NEW1.LT.1)G0  TO  341 

IF(RX.LT.RMIN.OR.RX.GT.RMAX)GO  TO  340 

IF(T.LT.TMIN.OR.T.GT.TMAX)GO  TO  340 
C  FOR  MORE  ACCURATE  CALCULATION  SET  THESE  MOLECULES  IN  OUTPUT  FLOWS 
C 
C  DEFINE  THE  VELOCITY  COMPONENTS 

CALL  RANDU(P) 

B=VTER1XSQRT(-AL0GCP)) 

CALL  RANDU(P) 

BB=2.XPI*P 

VELX=VELX+BXCOS(BB) 

VELY=VELY+B*SIN(BB) 

CALL  RANDU(P) 

VELZ=VTER1*SQRTC-AL0G(P)) 
C  DEFINE  THE  NEW  MOLECULE  TO  BE  ACTIVATED 

DO  325  IACT=1,NM0L1 

IF(P1C4,IACT).EQ.-99.)G0  TO  326 

325  CONTINUE 

C  IF  THERE  IS  NO  ROOM  FOR  ADDITIONAL  MOLECULE  PRINT  'ALARM' 
IFCIACT.GE.NMOLDGO  TO  3004 

326  P1(4,IACT)=RX 
P1(5,IACT)=T 
P1(1,IACT)=VELX 
P1(2,IACT)=VELY 
P1(3,IACT)=VELZ 

340  CONTINUE 

C* REPEAT 'PROCEDURE* FOR* SPECIES  2 

341  IF(NEW2.LT.1)G0  TO  351 
DO  350  12=1, NEW2 

CALL  RANDU(P) 

ATIME=P*DTM 

CALL  RANDU(P) 

RSTART=RMIN+(FLOAT(I-l)+P)XDR 

VELX=FWP2(2,I,KR) 

VELY=FWP2C3,I,KR) 

TEM=FWP2(4,I,KR) 
C  THERMAL  VELOCITY 

VTER2  =  SQRT(2.XB0LTZXTEM/SPECC2,D) 

XX=RSTARTXCOS(ALFA)+VELXXATIME 

YY=RSTARTXSIN(ALFA)+VELYXATIME 

RX=SQRT(XXXX2+YYXX2) 

T=ATAN(YY/XX) 

IF(RX.LT.RMIN.OR.RX.GT.RMAX)GO  TO  350 

IF(T.LT.TMIN.OR.T.GT.TMAX)GO  TO  350 
C  FOR  MORE  ACCURATE  RESULTS  SET  THESE  MOLECULES  IN  OUTPUT  FLOWS 

C  DEFINE  VELOCITY  COMPONENTS 

CALL  RANDU(P) 

B=VTER2XSQRT(-AL0G(P)) 

CALL  RANDU(P) 

BB=2.XPIXP 

VELX=VELX+BXCOS(BB) 

VELY=VELY+BXSIN(BB) 

CALL  RANDU(P) 

VELZ=VTER2XSQRT(-AL0G(P)) 
C  FIND  A  NEW  MOLECULE  TO  BE  ACTIVATED 

DO  345  IACT=1,NM0L2 

IF  (P2(4,IACT).EQ.-99.)G0  TO  346 

345  CONTINUE 

C  IF  THERE  IS  NO  PLACE  THEN  PRINT  ALARM 
IF(IACT.EQ.NM0L2)G0  TO  3004 

346  P2(4,IACT)=RX 
P2(5,IACT)=T 
P2(1,IACT)=VELX 
P2(2,IACT)=VELY 
P2(3,IACT)=VELZ 

350  CONTINUE 
C 
C  REPEAT  PROCEDURE  FOR  SPEC-3 

C 

C  INPUT  E  MOLECULES 
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351  IF  (TMIN.LE.O.)GO  TO  390  SIM06490 

C  SIM06500 

C  SIM06510 

ANEW1=F0W1C1,I,KR,KS+1)  SIM06520 

C  SIM06530 

C  SIM06540 

NEW1=ANEW1  SIM06550 

REM=ANEW1-NEW1  SIM06560 

CALL  RANDU(P)  5IM06570 

IF(P.LT.REM)NEW1=NEW1+1  SIM06580 

C  SIM06590 

ANEW2=F0W2C1,I,KR,KS+1)  SIM06600 

NEW2=ANEW2  SIM06610 

REM=ANEW2-NEW2  SIM06620 

CALL  RANDU(P)  SIM06630 

IFCP.LT.REM)NEW2=NEW2+1  SIM06640 

C  SIM06650 

C  SIM06660 

IF(NEW1.LE.O)GO  TO  370  SIM06670 

C  SIM06680 

DO  362  11=1, NEW1  SIM06690 

CALL  RANDU(P)  SIM06700 

T=TMIN+DDALFA*P  SIM06710 

CALL  RANDUCP)  SIM06720 

RNEW=RMIN+(FL0AT(I-1)+P)*DR  SIM067  30 

CALL  RANDU(P)  SIM06740 

TEM=F0W1(4,I,KR,KS+1)  SIM06750 

VTER1=SQRT(2.XB0LTZ*TEM/SPEC(1,1))  SIM0676  0 

B=VTER1*SQRT(-AL0G(P))  SIM06770 

CALL  RANDU(P)  SIM06780 

BB=2.*PI*P  SIM06790 

VELX=F0W1(2,I,KR,KS+1)+BXC0S(BB)  SIM068  0  0 

VELY=F0W1(3,I,KR,KS+1)+BXSINCBB)  SIM06810 

CALL  RANDU(P)  SIM06820 

VELZ=VTER1XSQRT(-AL0GCP))  SIM06830 

C  FIND  A  NEW  MOLECULE  TO  BE  ACTIVATED  SIM06840 

DO  364  IACT=1,NM0L1  SIM06850 

IF(P1(4,IACT).EQ.-99.)G0  TO  366  SIM06860 

364  CONTINUE  SIM06870 

C  IF  THERE  IS  NO  ROOM  FOR  ADDITIONAL  MOLECULE  THEN  PRINT  'ALARM1  SIM06880 

IFCIACT.EQ.NMOLDGO  TO  3004  SIM06890 

366  P1(4,IACT)=RNEW  SIM06900 

P1C5,IACT)=T  SIM06910 

P1C1,IACT)=VELX  SIM06920 

P1(2,IACT)=VELY  SIM06930 

362  P1C3,IACT)=VELZ  SIM06940 

370  CONTINUE  SIM06950 

C  SIM06960 

C  REPEAT  FOR  SPECIES  2  SIM06970 

IF(NEW2.LE.0)GO  TO  390  SIM06980 

DO  382  12=1, NEW2  SIM06990 

CALL  RANDU(P)  SIM07000 

T=TMIN+DDALFA*P  SIM07010 

CALL  RANDU(P)  SIM07020 

RNEW=RMIN+CFLOAT(I-l)+P)*DR  SIM07  030 

CALL  RANDU(P)  SIM07040 

TEM=F0W2(4,I,KR,KS+1)  SIM07050 

VTER2=SQRT(2.*B0LTZ*TEM/SPECC2,1))  SIM07  06  0 

B=VTER2*SQRT(-AL0G<P))  SIM07070 

CALL  RANDU(P)  SIM07080 

BB=2.*PIXP  SIM07090 

VELX=F0W2(2,I,KR,KS+1)+B*C0S(BB)  SIM07100 

VELY=F0W2(3,I,KR,KS+1)+B*C0S(BB)  SIM07110 

CALL  RANDUCP)  SIM07120 

VELZ=VTER2XSQRT(-AL0G(P))  SIM07130 

C  FIND  A  NEW  MOLECULE  TO  BE  ACTIVATED  SIM07140 

DO  384  IACT=1,NM0L2  SIM07150 

IF(P2(4,IACT) .EQ.-99)GO  TO  386  SIM07160 

384  CONTINUE  SIM07170 

C  IF  THERE  IS  NO  ROOM  FOR  ADDITIONAL  MOLECULE  THEN  PRINT  ALARM  SIM07180 

IFCIACT.EQ.NM0L2)  GO  TO  3004  SIM07190 

386  P2(4,IACT)=RNEW  SIM07200 


P2C5,IACT)=T 

P2C1,IACT)=VELX 

P2C2,IACT)=VELY 
382   P2C3,IACT)=VELZ 

C 

C  REPEAT  FOR  SPEC-3 
C 

390  continue 
c'input' 's* 'molecules' 

DO    460    J  =  1,NAD 

IF(KR.LT.2)G0    TO    460 

ANEW1=FNN1C1,J,KR-1) 

IFCANEW1.LE.  .0O00DGO    TO    420 

NEW1=ANEW1 

REM=ANEW1-NEW1 

CALL  RANDU(P) 

IFCP.LT.REM)NEW1=NEW1+1 

DO  402  11  =  1, NEW1 

CALL  RANDU(P) 

T=TMIN+DDALFAX(P+FL0AT(J-1)) 

CALL  RANDUCP) 

RNEW=RMIN+P*DR 

CALL  RANDU(P) 

TEM=FNN1C4, J,KR-1) 

VTER1  =  SQRT(2.*B0LTZ*TEM/SPEC(1,D) 

B=VTER1*SQRTC-AL0GCP)) 

CALL  RANDUCP) 

BB=2.*PI*P 

VELX=FNN1C2,J,KR-1)+BXC0S(BB) 

VELY=FNN1C3,J,KR-1)+B*SINCBB) 

CALL    RANDUCP) 

VELZ=VTER1*SQRTC-AL0GCP)) 
C    FIND   A    NEW   MOLECULE   TO    BE   ACTIVATED 

DO    404    IACT  =  1,NM0L1 

IFCP1C4,IACT).EQ.-99.)G0   TO    406 
404    CONTINUE 

IFCIACT.EQ.NMOLDGO   TO    3004 
406    P1C4,IACT)=RNEW 

P1C5,IACT)=T 

P1C1,IACT)=VELX 

P1C2,IACT)=VELY 
402  P1C3,IACT)=VELZ 
420  CONTINUE 
C  REPEAT  FOR  SPECIES  2 

ANEW2=FNN2C1,J,KR-1) 

IFCANEW2.LT.  .OOOODGO    TO    440 

NEW2=ANEW2 

REM=ANEW2-NEW2 

CALL  RANDUCP) 

IFCP.LT.REM)NEW2=NEW2+1 

DO  422  12=1, NEW2 

CALL  RANDUCP) 

T=TMIN+DDALFAXCP+FL0ATCJ-1)) 

CALL  RANDUCP) 

RNEW=RMIN+P*DR 

CALL  RANDUCP) 

TEM=FNN2C4,J,KR-1) 

VTER2=SQRTC2.*B0LTZXTEM/SPECC2,1)) 

B=VTER2*SQRTC-AL0GCP)) 

CALL  RANDUCP) 

BB=2.*PI*P 

VELX=FNN2C2,J,KR-1)+B*C0SCBB) 

VELY=FNN2C3,J,KR-1)+BXSINCBB) 

CALL  RANDUCP) 

VELZ=VTER2*SQRTC-AL0GCP)) 
C  FIND  A  NEW  MOLECULE  TO  BE  ACTIVATED 

DO  424  IACT=1,NM0L2 

IFCP2C4,IACT).EQ.-99)G0  TO  426 
424  CONTINUE 

IFCIACT.EQ.NM0L2)G0  TO  3004 
426  P2C4,IACT)=RNEW 
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P2C5,IACT)=T  SIM07930 

P2(1,IACT)=VELX  SIM07940 

P2(2,IACT)=VELY  SIM07950 

422  P2C3,IACT)=VELZ  SIM07960 

440  CONTINUE  SIM07970 

C SIM07  98  0 

C  REPEAT  FOR  SPEC-3  SIM07990 

C SIM08000 

C  INPUT  'N'  MOLECULES SIM08010 

C  SPEC-1  SIM08020 

ANEW1=FSN1(1,J,KR+1)  SIM08030 

IF(ANEW1.LE.  .0O00DGO  TO  450  SIM08040 

NEW1=ANEW1  SIM08050 

REM=ANEW1-NEW1  SIM08060 

CALL  RANDU(P)  SIM08070 

IF  (P.LT.REM)NEW1=NEW1+1  SIM08080 

DO  442  11=1, NEW1  SIM08090 

CALL  RANDU(P)  SIM08100 

T  =  TMIN  +DDALFA*(P+FL0AT(J-1))  SIM08110 

CALL  RANDU(P)  SIM08120 

RNEW=RMAX-P*DR  SIM08130 

CALL  RANDU(P)  SIM08140 

TEM=FSN1(4,J,KR+1)  SIM08150 

VTER1=SQRTC2.XB0LTZXTEM/SPEC(1,1))  SIM0816  0 

B=VTER1*SQRT(-AL0G(P))  SIM08170 

CALL  RANDU(P)  SIM08180 

BB=2.XPI*P  SIM08190 

VELX=FSN1(2,J,KR+1)+BXC0S(BB)  SIM08200 

VELY=FSN1(3,J,KR+1)+BXSIN(BB)  SIM08210 

CALL  RANDU(P)  SIM08220 

VELZ=VTER1XSQRT(-AL0GCP))  SIM08230 

C  FIND  A  NEW  MOLECULE  TO  BE  ACTIVATED  SIM08240 

DO  444  IACT=1,NM0L1  SIM08250 

IF(PK4,IACT).EQ.-99.  )G0  TO  446  SIM08260 

444  CONTINUE  SIM08270 

IFCIACT.EQ.NMOLDGO  TO  3004  SIM08280 

446  P1(4,IACT)=RNEW  SIM0S290 

P1(5,IACT)=T  SIM08300 

P1(1,IACT)=VELX  SIM08310 

P1(2,IACT)=VELY  SIM08320 

442  P1(3,IACT)=VELZ  SIM08330 

450  CONTINUE  SIM08340 

C SIM08350 

C  INPUT  N  MOLECULES  SPEC-2  SIM08360 

ANEW2=FSN2(1, J,KR+1)  SIM08370 

IFCANEW2.LE.  .OOOODGO  TO  460  SIM08380 

NEW2=ANEW2  SIM08390 

REM=ANEW2-NEW2  SIM08400 

CALL  RANDU(P)  SIM08410 

IF(P.LT.REM)NEW2=NEW2+1  SIM08420 

DO  452  12=1, NEW2  SIM08430 

CALL  RANDU(P)  SIM08440 

T=TMIN+DDALFA*(P+FL0ATCJ-1))  SIM08  450 

CALL  RANDU(P)  SIM08460 

RNEW=RMAX-PXDR  SIM08470 

CALL  RANDU(P)  SIM08480 

TEM=FSN2(4,J,KR+1)  SIM08490 

VTER2  =  SQRT(2.XB0LTZXTEM/SPEC(2,D)  SIM08  50  0 

B=VTER2*SQRT(-AL0G(P))  SIM08510 

CALL  RANDU(P)  SIM08520 

BB=2XPIXP  SIM08530 

VELX=FSN2(2,J,KR+1)+BXC0S(BB)  SIM08540 

VELY=FSN2(3,J,KR+1)+B*SIN(BB)  SIM08550 

CALL  RANDU(P)  SIM08560 

VELZ=VTER2*SQRT(-AL0G(P))  SIM08  57  0 

C  FIND  A  NEW  MOLECULE  TO  BE  ACTIVATED  SIM08580 

DO  454  IACT=1,NM0L2  SIM08590 

IF(P2(4,IACT) .EQ.-99)G0  TO  456  SIM08600 

454  CONTINUE  SIM08610 

IF(IACT.EQ.NM0L2)G0  TO  3004  SIM08620 

456  P2(4,IACT)=RNEW  SIM08630 

P2(5,IACT)=T  SIM08640 
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c 
c 

Cx 
Cx 

c 

5 
Cx 


P2(1,IACT)=VELX 
P2(2,IACT)=VELY 

452  P2(3,IACT)=VELZ 

460  CONTINUE 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

REARRANGE  MOLECULES  IN  THEIR  CELLS 
INITIALIZATION 

000  CONTINUE 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

KIP=NM0L1+NM0L2+NM0L3 
DO  1001  KADD=1,KIP 

001  IP(KADD)=0 
KADD=0 
N1A  =  0 

N2A  =  0 
N3A  =  0 
N1A,N2A,N3A  ARE  THE  NUMBER  OF  ACTIVE  MOLECULES  (COUNTED 

DO  1100  1=1, NRD 
DO  1100  J=1,NAD 

DO  1005  K  =  l,4 
005  ICCK,I,J)=0 

SET  SPECIES  1 
N1C  =  0 
N2C  =  0 
N3C  =  0 
NTC  =  0 

DO  1020  K1=1,NM0L1 
IF(P1(4,K1) .EQ.-99.)G0  TO  1020 
RLOC  =  ABS(CC19,I,  J)-P1(4,K1)) 
TLOC=ABS(CC20,I,J)-P1(5,K1)) 
IF(RL0C.GT.DRX.5)G0  TO  1020 
IF(TL0C.GT.DDALFAX.5)G0  TO  1020 

N1C=N1C+1 

NTC=NTC+1 

N1A=N1A+1 

IF(NTC.EQ.1)IC(4,I,J)=KADD 

IC(1,I,J)=IC(1,I,J)+1 

KADD=KADD+1 

IP(KADD)=K1 
.020  CONTINUE 
SET  SPEC-2  MOLECULES 

DO  1040  K2=1,NM0L2 

IFCP2(4,K2).EQ.-99.)G0  TO  1040 

RL0C=ABS(C(19,I,J)-P2(4,K2)) 

TL0C=ABS(C(20,I,J)-P2(5,K2)) 

IF(RL0C.GT.DRX.5)G0  TO  1040 

IF(TL0C.GT.DDALFAX.5)G0  TO  1040 


NEXT) 


N2A=N2A+1 

NTC=NTC+1 

IF(NTC.EQ.1)IC(4,I,J)=KADD 

IC(2,I,J)=ICC2,I,J)+1 

KADD=KADD+1 

IP(KADD)=K2 
1040  CONTINUE 
1100  CONTINUE 

WRITE(6,1021)N1A,N2A 
1021    FORMATC       NUMBER    OF   ACTIVE 


MOLECULES    SPEC1=    ',15, 


C 
Cx 

C 
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DO  999  1=1, NRD  SIM09370 

DO  999  J=1,NAD  SIM09380 

999  NCOL(I,J)=0  SIM09390 

C  SIM09400 

DO  900  1=1, NRD  SIM09410 

DO  900  J=1,NAD  SIM09420 

C  SIM09430 

NM(1)  =  IC(1,I,J)  SIM09440 

NM(2)=ICC2,I,J)  SIM09450 

NM(3)=IC(3,I,J)  SIM09460 

NT0T=NM(1)+NM(2)+NM(3)  SIM09470 

DO  900  L=l,3  SIM09480 

DO  900  M=l,3  SIM09490 

C  NUMBER  OF  SPECIES  IN  THE  PROGRAM  IS  2  (3)  SIM09500 

KV=(L-1)*3+M  'SIM09510 

KT=KV+9  SIM09520 

C  KV,KT  ARE  THE  ADDRESSES  OF  RELATIVE  VELOCITY  AND  COLLISION  TIMERS     SIM09530 

920  IF(C(KT,I,J).GE.TIME)  GO  TO  900  SIM09540 

C  C(KT,I,J)IS  THE  INTEGRATED  TIME  FOR  L-M  COLLISION  SIM09550 

KSEL=0  3IM09560 

KREJ=0  SIM09570 

IF(NM(L).GT.l.AND.NM(M).GT.l)GO  TO  912  SIM09580 

C  NO  COLLISIONS  ARE  CALCULATED  IF  THERE  ARE  NO  MOLECULES                SIM09590 

911  C(KT,I,J)=C(KT,I,J)+DTM  SIM09600 
GO  TO  900  SIM09610 

C  SELECT  NOW  THE  MOLECULES  FOR  COLLISION  SIM09620 

912  IF  (KSEL.GE.IOO)GO  TO  911  SIM09630 
CALL  RANDU(P)  SIM09640 
M0L1=PXNM(L)+. 999999  SIM09650 
IF(MOL1.EQ.O)MOL1=1  SIM09660 

C  SIM09670 

CALL  RANDU(P)  SIM09680 

M0L2=P*NM(M)+. 999999  SIM09690 

IF(MOL2.EQ.0)MOL2=l  SIM09700 

C  SIM09710 

KSEL=KSEL+1  SIM09720 

C  SIM09730 

C  CHECK  IF  THE  SAME  MOLECULE  HAS  BEEN  SELECTED  TWICE                    SIM09740 

IF(L.EQ.M.AND.M0L1.EQ.M0L2)G0  TO  912  SIM09750 

C SIM0976  0 

C  FIND  THE  ACTUAL  ADDRESSES  OF  THE  SELECTED  MOLECULES  SIM09770 

IF(L.EQ.1)K1=0  SIM09780 

IF(L.EQ.2)K1=NMC1)  SIM09790 

IF(L.EQ.3)K1=HM(1)+NM(2)  SIM098  0  0 

IF(M.EQ.1)K2=0  SIM09810 

IF(M.EQ.2)K2=NM(1)  SIM09820 

IF(M.EQ.3)K2=NM(1)+NM(2)  SIM09830 

KAD1=M0L1+K1+IC(4,I,J)  SIM09840 

KAD2=M0L2+K2+IC(4,I,J)  SIM09850 

C  KAD1,KAD2  ARE  THE  LOCATION  OF  SELECTED  MOLECULES  IN  IPC   )  SIM09860 

MAD1=IP(KAD1)  SIM09870 

MAD2=IP(KAD2)  SIM09880 

C  MAD1,MAD2  ARE  THE  ACTUAL  ADDRESSES  OF  THE  SELECTED  MOLECULES  (THE  SIM09890 

C  INDICATION  OF  WHAT  SPECIES  THEY  ARE  HAS  BEEN  DEFINED  HERE)             SIM09900 

DO  930  N=l,3  SIM09910 

IF(L.EQ.1)VN1=P1(N,MAD1)  SIM09920 

IF(L.EQ.2)VN1=P2(N,MAD1)  SIM0993  0 

C      IF(L.EQ.3)VN1=P3(N,MAD1)  SIM09940 

C  SIM09950 

IF(M.EQ.1)VN2=P1(N,MAD2)  SIM0996  0 

IF(M.EQ.2)VN2=P2(N,MAD2)  SIM0997  0 

C      IF(M.EQ.3)VN2=P3CN,MAD2)  SIM09980 

C  SIM09990 

930  VRC(N)=VN1-VN2  SIM10000 

C  VRC(3)  CONTAIN  THE  THREE  RELATIVE  VELOCITY  COMPONENTS                 SIM10010 

VR=SQRT(VRC(1)XVRC(1)+VRC(2)*VRC(2)+VRC(3)*VRC(3))  SIM10  020 

C  VR  IS  THE  RELATIVE  SPEED  IN  A  SPECIFIC  COLLISION  SIM10030 

IF(C(KV,I,J).LT.VR)C(KV,I,J)=VR  SIM10  040 

C  LAST  STATEMENT  RESETS  THE  MAXIMUM  RELATIVE  VELOCITY  FOR  FURTHER       SIM10050 

C  CALCULATIONS  SIM10060 

C  SIM10070 

IF(KREJ.GT.100)GO  TO  911  SIM10080 
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C  SIM1009I 

CALL    RANDU(P)  SIM1O10! 

AVR  =  VR/C(KV,I,J)  SIM1011I 

KREJ  =  KREJ  +  1  SIM1012I 

IFCAVR.LT.P)GO    TO    912  SIM1013I 

C    LAST    STATEMENT    REJECTS   THE   CALCULATED   COLLISION  SIM1014I 

C  SIM1015I 

C  NOW  A  SPECIES  L-M  COLLISION  HAS  BEEN  SELECTED  SIM1016J 

C  SIM10171 

C  CALCULATE  NOW  THE  PROBABILITY  THAT  SUCH  A  COLLISION  WILL  BE  COUNTED    SIM1018I 

C  FOR  THE  L  AND  M  SPECIES  RESPECTIVELY  SIM1019I 

LP  =  1  SIM1020( 

LM=1  SIM102K 

CALL  RANDUCP)  SIM1022I 

ANM=FLOATCNMCL))/FLOATCNMCM))  SIM1023( 

IFCANM.GT.l. )GO  TO  950  SIM10241 

IF(ANM.LT.P)MP  =  0  SIM1025C 

GO   TO    955  SIM1026C 

C  S1M1027C 

C  SIM1O280 

950    ANM=1./ANM  SIM1O290 

IF(ANM.LT.P)LP=0  SIM1O3O0 

Z  SIM1031C 

C  SIM1032C 

955   CXS  =  PI*CSPECCL,2)+SPECCM,2))X*2/<».  SIM1033C 

ALP  =  LP  SIM103^»C 

ALM  =  LM  SIM1035C 

C  SIM1036C 

V0LUME=REGC4,KR,KS)*V0LCI)  SIM1037C 

C CHECK  SIM1038C 

DNL=FLOAT(NM(L))/VOLUME  SIM1039C 

DNM=FLOAT(NM(M))/VOLUME  SIMIO^OC 

C    USE    EQ.10.3  SIM1041C 

TOCCL,M)=ALP/CCXS*DNMXVRXNMCL))+ALM/CCXSXDNL*VRxNMCM))  SIM1042I 

C  SIM1043( 

C    SET    THIS    VALUE    INTO    C(KT,I,J)  SIM1044( 

CCKT,I,J)=CCKT,I,J)+TOCCL,M)  SIM1045( 

C.  . .  SIM10461 

C    SAMPLE   THIS    COLLISION  SIM1047I 

NCOLCI,J)=NCOLCI,J)+l  SIM10481 

Z  SIM1049I 

C    FIND    RELATIVE   MASSES  SIM1050I 

CALL    RANDUCP)  SIM105K 

BB=1.-2.XP  SIM10521 

AA=SQRT(1.-BBXBB)  SIM1053I 

VRCC1)=BBXVR  SIM1054I 

CALL    RANDUCP)  SIM10551 

BB=2.XPIXP  SIM1056I 

VRCC2)=AA*C0SCBB)*VR  SIM1057I 

VRCC3)=AAXSINCBB)»VR  SIM1058I 

C  SIM1059I 

C  FIND  RELATIVE  MASSES  SIM1060I 

SM=SPECCL,1)+SPECCM,1)  SIM1061I 

RML=SPECCL,1)/SM  SIM1062I 

RMM=SPECCM,1)/SM  SIM1063I 

C  SIM1064I 

:xkxx»xxxxxxxxxxxxxxxxxxxxxxxxxxx«kxxxxxxxxx»xkxxxxxxxxxxxxxxxxxxxxkxxxxsIM1065I 

Z   CALCULATE  HERE  THE  ACTUAL  ADDRESSES  OF  COLLIDING  MOLECULES  AND  SET 

C  THEIR  NEW  VELOCITIES  CNEW  VELOCITY  COMPONENTS  ARE  ADDED  TO  THE 

C  VELOCITY  OF  THE  CENTER  OF  MASS  VCCM) 

DO  960  N  =  l,3 

1)V1=P1CN,MAD1) 
2)V1=P2CN,MAD1) 
3)V1=P3CN,MAD1) 
1)V2=P1CN,MAD2) 
2)V2=P2CN,MAD2) 
3)V2=P3CN,MAD2) 

VCCM=RMLXV1+RMMXV2 

VCCM1=VCCM+VRCCN)*RMM 

VCCM2=VCCM-VRCCN)XRML 
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IFCL.EQ 

IFCL.EQ 

C 

IFCL.EQ 

IFCM.EQ 

IFCM.EQ 

C 

IFCM.EQ 

C  CHANGE  IN  VELOCITY  IS  INPUTED  ONLY  IF  PROBABILITY  OF  COLLISION .GT . 1    SIM10810 

IFCLP.NE.DGO  TO  961  SIM10820 

IF(L.EQ.1)P1(N,MAD1)=VCCM1  SIM10830 

IFCL.EQ.2)P2(N,MAD1)=VCCM1  SIM10840 

C      IF(L.EQ.3)P3(N,MAD1)=VCCM1  SIM10850 

C  SIM10860 

961  IFCLM.NE.DGO  TO  960  SIM10870 

IF(M.EQ.1)P1(N,MAD2)=VCCM2  SIM1088  0 

IFCM.EQ.2)P2(N,MAD2)=VCCM2  SIM1089  0 

C      IF(M.EQ.3)P3(N,MAD2)=VCCM2  SIM10900 

C  SIM10910 

960  CONTINUE  SIM10920 

GO  TO  920  SIM10930 

900  CONTINUE  SIM10940 

WRITE(6,4545)T0CC1,1),T0CC1,2),T0CC2,1),T0C(2,2),TIME  SIM10950 

4545  FORMATC  TIME1 ,  5E10  .  3)  SIM10960 

C  SIM10970 

Cxxxx*x*x***x*xxx**END  OF  COLLISIONSxx*****xxxxxxxxxx*xxxxxxxxxxxxxxxxxxsiM1098  0 

C  NON  NEW  TEMPERATURES  MAY  BE  CALCULATED  IN  EACH  CELL  SrW10990 

C  AVERAGE  VELOCITY  IN  CELLS  SIM11000 

C  SIM11010 

WRITEC6,4547)  SIM11020 

4547  FORMATC «0     I     J   N1E    N2E      TEMPR1       TEMPR2        VAVX1  SIM11030 

1    VAVY1        VAVX2        VAVY2        NCOL',/)  SIM11040 

C  SIM11050 

DO  1110  1=1, NRD  SIM11060 

DO  1110  J=1,NAD  SIM11070 

N1E=IC(1,I,J)  SIM11080 

N2E=ICC2,I,J)  SIM11090 

KAD1=IC(4,I,J)  SIM11100 

KAD2=IC(4,I,J)+N1E  SIM11110 

IFCN1E.LT.DG0  TO  1112  SIM11120 

C  SIM11130 

VX1=0.  SIM11140 

VX2=0.  SIM11150 

VY1=0.  SIM11160 

VY2=0.  SIM11170 

C  SIM11180 

DO  1111  IM=1,N1E                               .  SIM11190 

MAD1=KAD1+IM  SIM11200 

IAD=IP(MAD1)  SIM11210 

VX1=VX1+P1(1,IAD)  SIM11220 

1111  VY1=VY1+P1(2,IAD)  •   SIM11230 

1112  IF(N2E.LE.0)GO  TO  1110  SIM11240 
C  SIM11250 

DO  1115  IM=1,N2E  SIM11260 

MAD2=MAD1+IM  SIM11270 

IAD=IP(MAD2)  SIM11280 

VX2=VX2+P2(1,IAD)  SIM11290 

1115  VY2=VY2+P2C2,IAD)  SIM11300 
C  SIM11310 
C  THE  AVERAGE  VELOCITIES  IN  THE  CELL  WILL  RESULT-  SIM11320 

VAVX1=VX1/N1E  SIM11330 

VAVY1=VY1/N1E  SIM11340 

VAVX2=VX2/N2E  SIM11350 

VAVY2=VY2/N2E  SIM11360 

IFd.NE.DGO  TO  1116  SIM11370 

C  SIM11380 

FSN1C2,J,KR)=VAVX1+FSN1(2,J,KR)  SIM1139  0 

FSN1C3, J,KR)=VAVY1+FSN1(3, J,KR)  SIM1140  0 

FSN2(2,J,KR)=VAVX2+FSN2(2,J,KR)  SIM11410 

FSN2(3,J,KR)=VAVY2+FSN2(3,J,KR)  SIM11420 

C  SIM11430 

1116  IF(I.NE.NRD)GO  TO  1117  SIM11440 
FNN1C2,  J,KR)=VAVX1+FNN1(2,  J,KR)  SIM11450 
FNN1(3,J,KR)=VAVY1+FNN1(3,J,KR)  SIM1146  0 
FNN2(2,J,KR)=VAVX2+FNN2(2,J,KR)  SI Ml 147  0 
FNN2C3,J,KR)=VAVY2+FNN2(3,J,KR)  SIM1148  0 

C  SIM11490 

1117  IFCJ.NE.DGO  TO  1118  SIM11500 
FENK2,I)=VAVX1  +  FEN1(2,I)  SIM11510 
FEN1(3,I)=VAVY1+FEN1(3,I)  SIM11520 
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FEN2(2,I)=VAVX2+FENZC2,I)  5IM11530 

FEN2C3,I)=VAVY2+FEN2(3,I)  SIM11540 

C  SIM11550 

1118  IFCJ.NE.NAD)GO  TO  1119  SIM11560 
FWN1C2,I)=VAVX1+FWN1(2,I)  SIM11570 
FWN1(3,I)=VAVY1+FWN1C3,I)  SIM11580 
FWN2(2,I)=VAVX2+FWN2(2,I)  SIM11590 
FWN2C3,I)=VAVY2+FWN2(3,I)  SIM11600 

1119  CONTINUE  SIM11610 
C  SIM11620 
C   THERMAL    VELOCITIES   AND   TEMPERATURES  SIM11630 

ENRG1=0.  SIM11640 

ENRG2=0.  SIM11650 

IF(N1E.LT.1)G0   T01131  5IM11660 

DO    1130    IM=1,N1E  SIM11670 

MAD1=KAD1+IM  SIM11680 

IAD=IP(MAD1)  SIM11690 

CX1=P1(1,IAD)-VAVX1  SIM11700 

CY1=P1C2,IAD)-VAVY1  SIM11710 

CZ1=P1(3,IAD)  S1M11720 

1130  ENRG1=ENRG1+CX1XCX1+CY1XCY1+CZ1*CZ1  SIM11730 

1131  IFCN2E.LT. DGO  TO  1150  SIM11740 
DO  1140  IM=1,N2E  SIM11750 
MAD2=KAD2+IM  SIM11760 
IAD=IP(MAD2)  SIM11770 
CX2=P2(1,IAD)-VAVX2  SIM11780 
CY2=P2C2,IAD)-VAVY2  SIM11790 
CZ2=P2(3,IAD)  SIM11800 

C  SIM11810 

C  SIM11820 

1140  ENRG2=ENRG2+CX2XCX2+CY2*CY2+CZ2*CZ2  SIM11830 

C  SIM11840 

1150  TEMPR1=(ENRG1/N1E)XSPEC(1,1)/(3.XB0LTZ)  SIM11850 
TEMPR2=(ENRG2/N2E)XSPEC(2,1)/C3.*B0LTZ)  SIM11860 

C  SIM11870 

IFCI.NE.DGO   TO    1151  SIM11880 

FSN1(4,J,KR)=TEMPR1+FSN1C4,J,KR)  SIM11890 

FSN2(4,J,KR)=TEMPR2+FSN2(4,J,KR)  SIM11900 

C  SIM11910 

1151  IF(I.NE.NRD)GO  TO  1152  SIM11920 
FNN1(4,J,KR)=TEMPR1+FNN1(4,J,KR)  SIM11930 
FNN2(4,J,KR)=TEMPR2+FNN2(4,J,KR)  SIM11940 

C  SIM11950 

1152  IFCJ.NE.DGO  TO  1153  SIM11960 
FWN1(4,I)=TEMPR1+FWN1(4,I)  SIM11970 
FWN2(4,I)=TEMPR2+FWN2(4,I)  SIM1198C 

C  SIM1199C 

1153  IF(J.NE.NAD)GO  TO  1154  SIM1200C 
FEN1(4,I)=TEMPR1+FEN1(4,I)  SIM1201C 
FEN2C4,I)=TEMPR2+FEN2C4,I)  SIM1202C 

1154  CONTINUE  SIM1203C 
C  SIM1204C 

WRITEC6, 4548)1, J, N1E, N2E, TEMPR1 , TEMPR2, VAVX1 , VAVY1 , VAVX2, VAVY2, NCOSIM1205C 

1L(I,J)  SIM1206C 

4548  FORMATC  • , 415, 6F12 . 3, 17 )  SIM1207C 

C  SIM1208D 

1110  CONTINUE  SIM1209C 

C  SIM1210C 

6000  CONTINUE  SIM1211C 

C  CALCULATE  AVERAGED  PARAMETERS,  WEIGHTED  BY  DFI  SIM1212C 

DO  6010  1=1, NRD  SIM1213C 

DO  6010  KPAR=1,4  SIM1214C 

F0W1(KPAR,I,KR,KS)=FWN1(KPAR,I)/(NIS)  SIM1215C 

F0W2(KPAR,I,KR,KS)=FWN2(KPAR,I)/(NIS)  SIM1216C 

F0E1(KPAR,I,KR,KS)=FEN1(KPAR,I)/(NIS)  SIM1217C 

6010  F0E2(KPAR,I,KR,KS)=FEN2(KPAR,I)/(NIS)  SIM1218C 

C  SIM1219C 

DO  6020  J=1,NAD  SIM1220C 

DO  6020  KPAR=1,4  SIM1221C 

FSNKKPAR,  J,KR)  =  FSN1(KPAR,  J,KR)/NIS  SIM12220 

FSN2CKPAR, J , KR) =FSN2( KPAR, J,KR)/NIS  SIM12230 

FNN1(KPAR,J,KR)=FNN1(KPAR,J,KR)/NIS  SIM12240 
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6020  FNN2(KPAR,J,KR)=FNN2(KPAR,J,KR)/NIS  SIM12250 

C  SIM12260 
C  TO  PREPARE  FLOWS  FOR  THE  NEXT  REGION  OR  SECTOR,  DIVIDE  FONO  BY  L0CALSIM12270 

C  DFI  AND  MULTIPLY  BY  NEXT  DFI  WHEN  STARTING  NEW  REGION  SIM12280 

C  SIM12290 

C  CALCULATE  HERE  THE  MEAN  FREE  PATH  AND  STORE  INTO  REG( , , )  SIM12300 

C  SIM12310 
C  STOP  PROGRAM  IF  FLOW  BECOMES  COLLISIONLESS  OR  NUMBER  DENSITY  IS  EQUAL  SIM12320 

C  TO  THE  AMBIENT  NUMBER  DENSITY.  STORE  F0E1,F0E2  IN  A  SEPARATE  FILE  TO   SIM12330 

C  BE  USED  IN  A  DIFFERENT  PROGRAM.  SIM12340 

C  SIM12350 

C  PRINT  AVERAGED  RESULTS  SIM12360 

WRITE(6,6029)  SIM12370 

6029  FORMATC  ',////)  SIM12380 
WRITE(6,6030)NIS,KR  SIM12390 

6030  FORMATC  AVERAGED  OUTPUT  FLOWS  AFTER1, 15,'  TIME  INCREMENTS  IN  REGISIM12400 
10NM5)  SIM12410 

WRITE(6,6031)  SIM12420 

6031  FORMATCO     I     FEN1(1,I)     FEN2(1,I)     FWN1(1,I)    FWN2( 1 , I) • SIM12430 
1)  SIM12440 

DO  6033  1=1, NRD  SIM12450 

WRITE(6,6  032)I,FEN1(1,I),FEN2(1,I),FWN1(1,I),FWN2C1,I)  SIM1246  0 

6032  FORMATC  « ,  15,  4F13  .  3)  SIM12470 

6033  CONTINUE  SIM12480 
WRITE  (6,6034)  SIM12490 

6034  FORMATCO     J     FNN1(1,J,KR)    FNN2(1,J,KR)    FSN1(1,J,KR)  FSNSIM12500 
12(1, J, KR)')  SIM12510 

DO  6036  J=1,NAD  SIM12520 

WRITE(6,6  037)J,FNN1(1,J,KR),FNN2(1,J,KR),FSN1(1,J,KR),FSN2(1,J,KR)SIM12530 

6037  FORMATC  ',I5,4F15.5)  SIM12540 

6036  CONTINUE  SIM12550 

C  SIM12560 

Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx  SI Ml 257  0 

C  SIM12580 

C  START  A  NEW  REGION  SIM12590 

C      IF(KR.EQ.IO)  GO  TO  7000  SIM12600 

C      KR=KR+1  SIM12610 

C      GO  TO  2000  SIM12620 

C  SIM12630 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM126  4  0 

C  SIM12650 

7000  CONTINUE  SIM12660 

C  SIM12670 

C  SIM12680 
Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx  SIM126  9  0 

C  SIM12700 

C  FIND  IF  FLOW  BECAME  COLLISIONLESS  IN  THE  WHOLE  SECTOR  SIM12710 

C  IF  POSITIVE,  STORE  F0E1,F0E2  IN  A  SEPARATE  FILE  AND  STOP  PROGRAM  SIM12720 

C  SIM12730 

C  PREPARE  DATA  FOR  THE  NEXT  SECTOR  SIM12740 

C      KR=1  SIM12750 

C      STOP  PROGRAM  IF  KS  WAS  BOUNDED  BY  THE  WALL  SIM12760 

C      KS=KS+1  SIM12770 

C      GO  TO  1000  SIM12780 

C  SIM12790 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM128  0  0 

C  SIM12810 
C  IF  THERE  IS  BACK  FLOW  (FWN1,FWN2)  NONZERO  CALCULATE  NEXT  ITERATION     SIM12820 

C      ITER=ITER+1  SIM12830 

C      KS=KR=1  SIM12840 

C      GO  TO  3000  SIM12850 

C  SIM12860 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM1287  0 

C  SIM12880 

WRITE(6,6099)  SIM12890 

6099  FORMATC1  DATA  FOR  TEN  MOLECULES  SPEC.2')  SIM12900 

DO  3003  1=1,10  SIM12910 

WRITE(6,3002)  P2( 1 , 1) , P2( 2, 1 ) , P2( 3, 1 )  , P2( 4, 1 ) , P2(5, 1 )  SIM12920 

3002  FORMATC  > , 4F10 . 3, E18 . 10 )  SIM12930 

3003  CONTINUE  SIM12940 
GO  TO  3009  SIM12950 

3004  WRITE(6,3005)  SIM12960 
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3009  CONTINUE  N0  PLACE  F0R  ADDITI0NAL  MOLECULES')  SIM12970 

cTnp  SIM12980 

FND  SIM12990 

c        U  SIM13000 

C  SIM13010 

SUBROUTINE  RANDU(P)  UmJ^?! 

COMMON  IX  I™?*??! 

IY  =  IXX65539  t?m  ™*2 

IF  (IY)  5,6,6  IrSilnlS 

5  IY_*  IY+2147«36«7  +  1  limllll 

P  =_P*  «656613E-9  2IM13080 

RETURN  liHfiJf! 

Fwn  SIM13110 

SIM13I20 
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B.5   Program  SIMUL  -  User's  Guide 

Preparation:   Run  AXSYI1  program.   From  its  output  data  evaluate: 
ALFA  -  The  averaged  angle  for  the  continuum  breakdown 

(P  -  0.05) 
TETA  -  Flow  direction  along  the  breakdown  limit. 

Flow  parameters  along  this  line  -  Pressure,  temperature,  velocity,  mean 
free  path. 

Input  Data: 

Radius  of  nozzle  ring  Rl 

Radial  size  of  a  cell  DR 

Maximum  radius  in  simulation  RP 

Angle  of  breakdown  limit  ALFA 

Flow  direction  along  (ALFA)  TETA 

Averaged  flow  velocity  along  (ALFA)  Vo 

Molecular  weight  of  each  species  Spec(I,l) 

Molecular  diameter  of  each  species  Spec  (1,2) 

Mean  free  path  along  (ALFA)  FPM 

Averaged  Temperature  (ALFA)  TEMP 

Time  increment  DTM 

Number  of  time  increments  NIS 

Options 

a.   Geometry 

The  program  is  designed  to  run  for  axisymmetric  ring  flow. 

For  a  two  dimensional  planar  flow  the  molecular  motion,  collisions  and 

flow  calculation  remain  unchanged.   The  flow  cross  section  remain 
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unchanged.   Instead  of  the  angle  DFI  the  two  dimensional  flow  requires 
the  definition  of  the  width  of  each  region  (or  cell).   To  keep  the 
number  of  molecules  within  reasonable  computational  limits  this  size 
has  to  decrease  in  the  same  manner  as  DFI. 

The  part  of  the  program  which  has  to  be  changed  for  this  purpose  is 
lines  230  to  250. 
b.   Wall  flux  calculations: 

The  program  may  run  for  the  whole  molecular  region  resulting  the  flux 
towards  the  wall  (this  part  needs  additional  debugging).   However  in 
order  to  make  it  more  efficient  we  may  stop  the  program  at  a  sector 
where  the  flow  becomes  collisionless .   The  remaining  part  of  the  flow 
may  be  regarded  either  as  collisionless  or  if  we  define  a  much  larger 
size  of  cells  we  may  calculate  the  molecular  collisions  on  this 
basis . 
This  part  need  additional  analysis  and  programming. 

Execution  Commands 

Without  additional  changes  the  program  runs  under  WATFIV  compiler  using 

the  following  command: 

WATFI  AXSYM  *  (XTYPE) 

Further  developments  will  be  required  to  run  the  program  for  each  sector 

separately  and  the  output  intermediate  results  on  the  mass  storage. 
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B.6     The  Influence  of  the  Ambient  Gas 

The  temperature,  pressure  and  density  of  the  ambient  gas  are  shown  in 
Table  1.   Because  its  temperature  is  much  higher  than  in  the  jet  gas,  the 
thermal  velocity  of  ambient  gas  molecules  is  much  higher  than  jet  molecules, 
the  following  contribution  may  be  expected  due  to  the  ambient  gas: 

a.  Collisions  between  the  "hot"  ambient  molecules  with  the  "cold"  jet 
molecules  may  cause  an  increase  in  the  dissipation  in  the  outer 
layer  of  the  jet  and  increase  in  the  flux  towards  the  walls. 

b.  For  higher  ambient  pressures,  those  collisions  become  rare  because 
of  the  low  number  density  therefore,  the  influence  of  the 
collisions  may  decrease. 

The  only  way  to  evaluate  the  influence  of  these  two  controversial  factors 
is  by  an  additional  simulation  program  to  be  designed  for  this  region.   The 
following  are  the  main  factors  to  be  included  in  this  program: 

Boundary  Conditions: 

a.  The  jet  side:   F0E1  and  F0E2  obtained  from  the  last  simulated  sector 
(SI1IUL)  supply  the  number  flux,  flow  velocity  componetns  and 
temperatures.   These  parameters  are  given  for  all  points  along  the 
radius  R(  I)  at  constant  distances  DR.   In  the  low  density  domain  the 
resolution  OR  is  much  too  large  compared  with  the  expected  aiean  free 
path.   A  different  mesh  has  to  be  designed  for  this  purpose.   It  is 
possible  that  one  cell  may  be  sufficient. 

b.  Far  Field  Condition:   The  boundary  conditions  where  the  gas  may  be 
regarded  "undisturbed"  by  the  jet  are  as  follows 
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-  Jet  gas  molecules  are  allowed  to  go  out  the  simulation  region 
(these  molecules  will  be  regarded  as  "lost"  molecules). 

-  Ambient  gas  is  allowed  to  enter  the  simulated  region  according 
with  their  thermal  velocity  and  number  density. 

c.   Solid  Wall  Bondaries .   The  solid  wall  may  be  assumed  to  have  a 

constant  temperature  Tw.  Incident  molecules  of  either  species  are 
reflected  back  from  the  wall.  Different  models  of  collisions  with 
the  wall  may  be  employed. 

-  Elastic  collisions  'specular  reglection'  (this  calculation  has  been 
included  in  SIMUL  program) 

-  Collisions  with  ideal  heat  transfer  (diffuse  reflections) 

-  Other  models  depending  on  the  materials  and  surface  parameters. 
The  collision  with  the  wall  was  included  as  an  internal  routine  in  the 

program  SIMUL.   If  the  general  molecular  simulation  contains  the  program 
proposed  here  the  collision  with  the  wall  will  be  omitted  from  SIMUL  and 
become  the  core  of  the  additional  collisionless  program.   Figure  25  shows  the 
low  density  (collisionless)  region  and  its  boundaries. 
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Figure   25.      The   Low  Density   Region   in    the   Jet 
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SUMMARY  OF  REPORT 

Algorithms  for  the  continuum  regime  and  for  the  region  where  molecular 
collisions  are  significant  have  been  developed. 

Program  AXSYM  contains  the  calculation  of  planar  jet  flow  and 
axisymmetric  ring  jet  flow.   This  program  supplies  data  for  the  limits  where 
the  continuum  approach  become  invalid  and  molecular  approach  should  be 
employed. 

Program  SIMUL  contains  the  molecular  simulation  for  the  axisymmetric  ring 
flow.   This  program  may  run  for  the  whole  molecular  region  to  result  in  the 
calculation  of  flux  towards  the  solid  wall.   For  a  more  efficient  simulation 
it  is  proposed  to  design  an  additional  program  for  the  collisionless  region 
where  ambient  gas  may  be  included. 

For  the  two  dimensional  flow,  program  SIMUL  may  be  used  after  changing 
the  definition  of  the  cell  geometry. 

To  run  the  whole  program  it  may  be  required  to  make  separate  runs  for 
each  sector  and  store  results  on  the  mass  storage. 
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